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Abstract 


This report concerns the prediction of the elastic moduli and the internal stresses within the 
unit cell of a fabric reinforced composite. In the proposed analysis no restrictions or assumptions 
are necessary concerning yam or tow cross-sectional shapes or paths through the unit cell but the 
unit cell itself must be a right hexagonal parallelepiped. All the unit cell dimensions are assumed 
to be gmall with respect to the thickness of the composite structure that it models.. 

The finite element analysis of a unit cell is usually complicated by the mesh generation 
problems and the non-standard, adjacent-cell, boundary conditions. This analysis avoids these 
problems through the use of preprogrammed boundary conditions and replacement materials (or 
elements). With replacement elements it is not necessary to match all the constituent material 
interfaces with finite element boundaries. Simple brick-shaped elements can be used to model the 
unit cell stmcture. The analysis predicts the elastic constants and the average stresses within each 
constituent material of each brick element. The application and results of this analysis are 
demonstrated through several example problems which include a number of composite 
microstructures. 
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I. Introduction 


A unit cell of fabric reinforced composite is any small, closed, polygonal volume of 
inhomogeneous material (often brick shaped) which, when reproduced and similarly aUgned, can 
be stacked , (side by side, top to bottom, and end to end) and joined together (as in solid brick 
construction) to approximate a variety of simple structural components whose minimum external 
dimensions are much larger than any unit cell dimension. Furthermore, it is desired that the 
thermo-mechanical response of the component and the unit ceU assembly be similar. A variety of 
different unit cells and analyses have successfully predicted fabric reinforced composite moduli 
(Ref. 1) and average thermal properties (Ref. 2) but the resolution of the detailed internal stress 
distribution within a unit cell has been more difficult. 

The ability to resolve the stresses within the unit cell of a fabric reinforced composite has 
at least three areas of applicability. Problems of crack growth within the microstructure are the 
most challenging of the three. The capacity to resolve the stress details must be very high in this 
appUcation. Another level of usefulness is the prediction of the initiation and propagation of 
yielding or plastic flow (usuaUy in the matrix phase) within the microstructure. This still requires 
a detailed knowledge of the internal stresses but it is not necessary to superimpose crack induced 
stresses on top of an already complicated stress field. A third, and much less demanding, level of 
usefulness is in material ranking and trade off studies. This level of engineering rates the 
likelihood of different fabric microstructures to perform satisfactorily in specific applications. 
Here the performance criteria can be quite simple and the demand for stress accuracy and detail 
can be significantly less than in the two prior applications. The large number of material and 
microstructural parameters available to the designer (or selector) of a fabric reinforced composite, 
coupled with the expense of experimentally characterizing these materials, makes initial screening 
by mechanical analysis more attractive. It was this application that was of most concern in the 

development of this analysis method. Numerical accuracy was clearly sacrificed to reduce 
modeling complexity in a manner consistent with material screening and comparison study 
requirements. 

The three-dimensional stresses within a unit cell of a fabric reinforced composite can be 
predicted by the application of a general purpose finite element code. However, the associated 
boundary conditions on the unit cell surface and the mesh generation problems can be difficult. 
The program described in this report avoids these difficulties through the use of preprogrammed 
boundary conditions and replacement elements. With replacement elements it is not necessary to 
match all the internal material interfaces with finite element boundaries. Thus, simple, uniform, 
parallelepiped elements can be applied to a unit cell structure whose boundaries are themselves a 



parallelepiped. Most of the common reinforcing microgeometries can be modeled with this shape 
of unit cell. The analysis predicts both the stresses (and strains) within each homogeneous 
element, and the average stress (and strain) within each dissimilar material contained in each 
replacement element. Conventional yield or failure criteria can then be applied to each material in 
each element, as in conventional stress analysis. 

The proposed analysis places no restrictions on fabric microgeometry within the unit cell 
except that the fibers all be continuous, the fiber packing within any tow remain relatively 
constant, and the microgeometry be deterministic. 

The key to the usefulness of this analysis is the performance of the replacement elements. 
This performance will be investigated for several sample problems of increasing complexity. 

These sample problems also help to explain the analysis and its application. The discussion begins 
with a simple one-dimensional tension bar problem. At this level the analysis seems almost trivial. 
The extension to two and three-dimensional problems is not trivial. In some of the sample 
problems the exact solution for the internal stresses is known. The plain weave unit cell is the most 
complex of the sample problems. For comparison, another numerical solution to this problem is 
available from an earlier study. 

The two and three-dimensional problems require a computer analysis. The final version 
of this numerical analysis, as it evolved from a sequence of programs directed at each sample 
problem, is a Fortran program written for the Sun Spark station I"*". All of the equations and 
derivations for the two and three-dimensional analyses, along with the program listing and 
input/output descriptions, appear in the Appendices. 

This analysis method and the related Fortran program, REPLACE, are considered to be an 
update of the earlier analysis program, FABNEW, which was developed about four years ago 
(Ref. 1). However, the earlier program has a thermal expansion prediction capability that could 
not be incorporated into REPLACE due to time and schedule limitations. 
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II . One Dimensional Analysis 

In this section, the appUcation and characteristics of replacement finite elements will be 

introduced at the simplest level, nainely one-dimensional ela^^ analysis. Through the example of 

a tension bar, the convergence of various finite element models for the elastic deformations will be 
investigated find compared to the known solution. The proposed replacement element analysis is 
also capable of predicting average stresses in each constituent material within each element The 
accuracy of these stress predictions are considered. There is no direct computational advantage to 
the use of replacement elements to model such a simple problem but it is instructive to initially 
consider the use of these elements at this elementary level. 


Sample Problem #1 

Consider the tension bar of Figure 1 in which the left hand half is made from a 
homogeneous isotropic material with modulus E and cross-sectional area A. The other half has 
the same cross-sectional area but the material is five times stiffer. From elcmentry considerations 
the total elongation of the bar (§ ) is given by the sum of the elongations of the two halves. 

where P is the axial load and L the total length of the bar. The axial stress ( 0 ^ ) and strain (£ ) in 
each material are given by 



where subscripts r, l designate right and left. 

The same results could also have been obtained using finite element analysis as long as one 
of the finite element nodes coincided with the material discontinuity. In that case all of the 
elements would be homogeneous and their stiffness matrices precise, as loi^ as the assumed 
displacement mode shapes included a constant and a linear term. The stiffness matrix [k] relevant 
to the axial forces and displacements at the end points of the bar is given by 


tk3 


5AE I “I 
3L -I , 
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If the material discontinuity does not coincide with a node point then one element will be 
inhomogeneous, as shown in Figure 1 , and the finite elonent solution will be an i^proximate one, 
as long as the assumed displacements are simple polynomials. The accuracy and convergence 
depends on the choice of mode shapes. For example, consider a linearly varying displacement 
within each element and an internal node placement at the 1/3 and 2/3 points along the bar length, 
as shown in Figure 1 . Each subsequent refinement of the finite element grid divides each prior 
element into three equal segments. The middle element of the model will always be 
inhomogeneous as the element size decreases. The stiffness matrices for the homogeneous elements 
arc given by 



where 1 is element length. 

The stiffness matrix for the single inhomogeneous element could be obtained from the 
general energy formula (Ref. 3) 

C kO = jjj B^D B dv = aTb^OBcIx (I) 

VOL 0 

where B is the strain/displacement matrix and dx (dv) is an increment of length (volume) along the 
bar. D is the local material stress/strain relation. Supercript t designates transpose of a matrix. 
The resulting inhomogeneous bar stiffness matrix is given by 



Figure 2 is a plot of the error in the bar elongation prediction as element size diminishes. The 
predicted end displacement approaches the known solution monitonically as the influence of the 
single inhomogeneous element error diminishes with element length. However, the error in the 
average strain of the center element persists at a high level (80%). This error can be reduced by 
resorting to higher order elements; but there is no accepted method for obtaining either the average 
or the detailed strains or stresses in the constituent materials within the inhomogeneous element. 
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Now consider a different approach to the same problem. Instead of q^plying the energy 
formula for the stiffness matrix, replace the inhomogeneous material with a fictitious homogeneous 
material that matches the axial response of the inhomogeneous materials, m center element is 
obviously a case of "stiffhess in series", for which an equivalent modulus ® can be obtained from 
the rule of mixtures for stiffnesses in series (Ref. 4): 


where V i stands for fractional length of the i th segment of the element . For the particular 
example at hand where E = 



If this equivalent modulus is used for the center element the exact solution results. What is better, 
the "stiffhess in series" model can be used to compute the average and local stresses and strains in 
the various materials of the inhomogeneous element from the nodal displacement solution. In 
particular, from Figure 1 , 



which is the correct result. 

This process of substituting equivalent homogeneous elements in place of inhomogeneous 
ones is termed the "replacement element" method. 

Of course, if the "element in series" results were known, a priori, there would have been no 
need to resort to a finite element solution. However, in more complicated two and three- 
dimensional problems, knowing the local solutions for series and parallel stiffness models is not 
equivalent to solving a global problem that involves their use in place of inhomogeneous elements. 
For example, if the tension bar of Figure 1 were part of a redundant truss problem a truss analysis 

would still be required. 

Tlie error inherent in the use of the general energy formula, in combination with a low 
order displacement mode shape assumption, arises from the formula's inability to distinguish 
between series and parallel stiffnesses. For one-dimensional problems, with linear displacement 
assumptions, the energy formula presumes a "stiffnesses in parallel" situation, whether that is the 

case or not. The introduction of higher order displacement modes permit the general energy 

formula to make the necessary distinction. However, for polynomial mode shapes and 



discontinuous material properties, the convergence rate improvements arc slow and detailed stress 
and strain determination problems remain. 
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III. Two Dimensional Analysis 

This section applies the concept of substituting replacement homogeneous elements 
in place of inhomogeneous ones at the generalized plane strain level of two-dimensional 
analysis. As in the one-dimensional case, the approach is flrst illustrated through a speciHc 
example for which the exact solution is easily obtained. 

There are several different notions that should be introduced in the transition from 
one to two dimensions. The first is the unit cell concept. Much of the earlier work (Ref. 5) 
on the resolution of detailed stress Helds in unidirectional materials (and laminates built up 
from unidirectional plies) used this type of idealization to make a large random 
microgeometry amenable to deterministic analysis. The unit cell approach looks for the 
simplest essential volume of composite microstructure from an analysis viewpoint. In 
two-dimensional analysis this selection is usually easy. Ref. 5 considered some convenient 
unit cells for square and hexagonally packed unidirectional composites. Each of the three 
sample problems in this section will begin by defining one or more unit cells for subsequent 
analysis. There are an infinite number of possible unit cells for a typical composite 
microstructure so the final choice is often somewhat personalized. The smallest unit cell is 
not always the most convenient one if the boundaries are non-rectangular. 

Another basic difference between one and two-dimensional problems is the 
mathematical nature of the replacement element idealization. In one dimension the material 
interfaces are discrete points. Continuity of normal stress and the geometric relationship 
between average element normal strain and average constituent normal strains are the only 
relevant concerns. In two-dimensional analysis the constituent material interfaces are 
assumed to be linear (or planar) with several local stress and strain components of concern. 

The physical nature of the replacement element process also changes from series 
and parallel bar or rod models to parallel plate models. The use of the general energy 
formula from Eq. 1 (as applied to a two-dimensional finite element) in combination with 
low order displacement mode shapes lead to the tacit assumption that each constituent 
material is arranged in a stacking of thin plates parallel to the plane of the analysis. The 
di ss i milar material plates have their thicknesses in proportion to their respective volume 
fractions in the element. In reality, the constituent material interfaces are not parallel to the 
analysis plane but normal to it. The replacement element process corrects this 
inconsistency by rotating the same stacking of plates 90^ about the material interface such 
that the tinal set of interfacial planes, between the parallel plates, preserves the original 
angle of the interface in the plane of the analysis. This procedure can only be applied to 
two constituent materials at a time whose interface is a single straight line in the 
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plane of the analysis. Thus, while the energy formula preserves only the constituent material 
volume fraction, the replacement element process preserves both the constituent volume fraction 
and the direction of the interface. Only the order or sequence of constituent material positioning 
across an interface is lost in the idealization. This process is best understood by considering the 
specific examples that follow. 


Sample Problem #2 

Figure 3 shows a laminated composite consisting of parallel bonded sheets of two 
different homogeneous isotropic materials. On a gross scale this assemblage of plates may be 
considered to be a composite material with a plane of isotropy parallel to the material interfaces. 
The principal axes of the composite are any pair of axes in the plane of isotropy with a third axis 
normal to that plane. In the principal axes, or natural coordinates of the composite, the elastic 
constants can be established from the application of elementry mechanics principals to the unit cell 
structure. Also, the same elementry model can be used to obtain the equations for the internal 
stresses in each contitutent material corresponding to any remotely applied state of uniform 
composite stress or strain. The elastic constants and the detailed stresses and strains can then be 
transformed into any global reference system; in particular, the one ^wn in Figure 3 where one of 
the natural coordinates correspond to the z-axis of the global reference system. 

The isotropic properties of the two sets of parallel plates can be chosen to match the 
properties of aluminum and epoxy from Table 1 . The volume fractions of both constituents are 
0.5. From elementry mechanics considerations the elastic constants of the composite, in the 
principal axes, can be obtained as follows. Consider the unit cell of Fig. 3 in the 1 ,2,3 coordinate 
system. From equilibrium and resolution of forces the average composite stresses ( (Jj , T; j ) are 
related to the constituent stresses ( CJ;', ) by 


a, = 
^2 = 
^3 = 

7-, 2 = 


a,*'VAi+crf'’VEp 

a/' = at 

at' Va1 + ‘^3'’Vep 
^Ai _ 

M2 “ M2 


where vj designates volume fraction of the i th constituent and Ep and Al designate epoxy and 
aluminum respectively. The corresponding strains ( € •, ^ Yij , , 7ij ) 

geometry and compatability as follows 
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t| - C| - t,| 

62 = +6 I** Vep 

^ _ ^Al . £ Ep 

7,2=T.2Vai + T,P Vep. 

These 12 equations plus the individual stress/strain laws for the two constituent materials form a 
system of 20 equations that can be solved for the composite stress/strain relation and the individual 
constituent stresses and strains corresponding to any applied composite stresses or strains (see 
Appendix A). From the composite stress/strain relations the composite elastic constants are 


El = Ej = 5.25 X 10<» psi 
E 2 » 1 . 39 x 106 psi 
Gi 2 = 0.354 X 10^ psi 
Vi2*V32*“ 0-325 
Vi 3 = 0.255 


From these principal values the engineering constants in another coordinate system, obtained by a 
rotation about the 3-axis of Figure 3, can be calculated from the appropriate 2-D transformation 
equations (given in Appendix B). In particular, for a rotation of 45® about the 3-axis of Figure 3 
the elastic constant are 

Ex = Ey»l.llxl0<>psi 

Ej = 5.25 X 10<» psi 
Gxy = 0.968 X 10^ psi 
Vxy“ 0.566 
Vx2 = Vyx = 0.067 

“nxyx * ^xy,y = -0.296 
%y,z'0034 


9 



For an average composite tensile stress of one psi in the x-direction (with all the other 
average stress components equal to zero) the stresses in the constituent materials arc given by 



Aluminum 

Epoxy 

Q 

1.103 

0.897 


0.103 

-0.103 


-0.251 

0.251 


0.103 

-0.103 


These stress and moduli predictions from elementry analysis are exact because they can be shown 
to satisfy all the local and global conditions of equilibrium and compatibility. 

As in the one-dimensional example, these results can also be obtained by conventional 
finite element analysis using various types of elements and grids. The unit cell can be analyzed in 
the principal coordinates of the material, as shown in Figure 4 , using rectangular or constant 
strain triangular elements without violating element material homogeniety. The applied unit stress 
in the x-dircction can be resolved into its components in the 1,2,3 coordinates of Figure 4 by either 
a Mohr's circle or the use of the stress transformation equations of Appendix B. The resulting 
composite moduli and constituent stress predictions can then be transformed back into the global 
x,y,z coordinate system. These results agree preeisely with the results of the elementry analysis. 

Alternatively, using the unit cell and grid of Figure 5 A, with constant strain triangular 
elements, the exact results can be obtained from homogeneous elements without the necessity of 
transforming the input and output from one coordinate system to another. 

It is interesting to also consider the application of inhomogeneous Hnite elements to the 
analysis of the same unit cell. Figure 5B shows this unit cell of the composite and one possible 
subdivision of the unit cell into rectangular elements. Some of the elements are homogeneous and 
some inhomogeneous. Using 4-node, isoparametric, brick elements (Ref 3); generalized plane 
strain analysis; the 25-node Hnite element grid shown in Figure 5B; and die general energy 
formula (Eq. 1) for the stiffness matrix of the inhomogeneous elements, the analysis 
overestimates the x and y moduli by almost 1 00%. Refinement of the grid leads to the moduli 
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estimates of Figure 6. The convergence is slow. Furthermore, there u no effective method of 
obtaining constituent material stresses within the inhoiiK^eneous dements. 

Now consider replacing the inhomogeneous elonents in this exanq)le problem with 
replacement elements. To make this substitution in two dimensions first consider a subelement of 
the inhomogeneous material, shown in Figure 7. The sides of tiiis subelement are dther parallel or 
normal to the material boundary plane. The volume fractions of the two materials are the same in 
the subelement as in the element that contains it. Assume that die replacement homogeneous 
material for the subelement and the whole element are the same. The derivation of Appendix A 
then can be applied to establish both the replacement homogeneous material moduli and the 
average constituent material stresses, once the average element strains are established. The physical 
nature of the homogeneous-inhomeogeneous replacement process is now evident. The 
inhomogeneous element of Figure 7 is replaced by a homogeneous composite element consisting of 
parallel plates bonded together in the same volume fraction as the inhomogeneous element and 
having the same orientation of the material interfaces. With the 25-node finite element grid the 
substitution is of the nature shown in Figure 8. For simplicity let the rectangular element 
stiffness matrix be made up of the sum of two constant strain triangular elements. (There is no 
n eed for higher order elements in this example.) The same replacement material substitution is 
done for both of the constant strain triangles that make up the rectangular element. The stress 
predictions for the constituent materials in the rectangular element are the average values from the 
two triangles. 

The results from the 25-node fmite element analysis are not the same as the exact solution 
for either the moduli or the constituent stresses. The Young's modulus in the loading direction is 
3 1% high as a result of the use of the replacement elements. This is a considerable improvement 
over the 100% error using the same finite element grid with the general energy formula for element 
stiffness. This error diminishes to less than 14% if the rectangular grid is changed from 4x4 to 8x8 
as shown in Figure 9. Since the replacement elanent analysis also provides constituent stresses it 
is of interest to compare the stresses in the 4x4 replacement elements to the known results. The 
following table makes this comparison. 




Phase 

Epoxy Phase 1 


Replacement 

Element 

Exact Result 

Replacement 

Element 

Exact Result 

^x(psi) 

0.906 


0.677 

0.898 

5y(psi) 

0.168 


-0.081 

-0.103 

^z(psi) 

-0.264 

-0251 

0.189 

0.251 

Txy(psi) 

0.155 

0.103 

-.075 

-0.103 


11 












The peak stresses from the replacement elements are about 20% lower than the exact values. 
Unfortunately, as in the one-dimensional case, these constituent stress errors do not diminish widi 
grid refinemoit. These errors must be reduced by the use of improved elements. The stresses in the 
homogeneous elements away from the replacement elements do converge rapidly to the exact 
results with increasing grid refinement. 

As was true in the tension bar example, the use of the general energy expression for the 
inhomogeneous element stiffness matrix, in combination with low order di^lacement mode 
shapes, favors an "elements in parallel" model of behavior rather than an "elements in series" 
model as is sometimes more appropriate. Figure 10 illustrates this tendency of an inhomogeneous 
plane stress element (by reference to a lattice or framework model). If the upper and lower halves 
of the element, as shown in Figure lOA , were made of dissimilar isotropic materials then good 
engineering judgment would dictate the lattice representation of Figure lOB, where lattice members 
that cross the material boundary are modeled as "elements in series" while those that do not cross 
the material boundary are simply homogeneous. The low order energy formula leads to a lattice 
structure of the type shown in Figure IOC. If there is not much difference between the stiffness of 
the constituent materials the two lattice models do not differ significantly. But if the constituents 
are very different, elastically, then the two models differ widely. 


Sample Problem #3 

This sample problem involves the determination of the extensional moduli and fiber/matrix 
stress concentrations for a unidirectional composite consisting of a square packed array of glass 
fibers in an epoxy matrix. These stiffnesses and stress concentrations are well establi^ied from 
several earlier micromechanics investigations. It will be shown that finite element analysis based 
on the substitution of orthotropic replacement elements for the inhomogeneous elements can yield 
q>proximately the same results for both moduli predictions and stress analysis even though the 
stresses within any constituent material in the unit cell model are not uniform. 

The specific problem concerns a 50% flber volume fraction of unidirectional E glass in an 
epoxy matrix. Figure 1 1 shows the square packed array of fiber cross-sections and a single unit 
cell of the composite. At most, only one quadrant of the unit cell needs to be analyzed due to 
structural and load symmetry. The constituent material properties are given in Table 1 . 

The 5x5 rectangular finite element grid of Figure 12 is superposed on the fiber/matrix 
geometry. The rectangular, generalized plane-strain, element stiffness matrices are formed from a 
pair of constant strain triangular elements, using the same replacement material properties in each 
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triangle of the rectangle. This leads to the material model of Figure 13 in \riiich the plate thickness 
and ^cing within each originally inhomogeneous element reflects the true constituent volume 
fractions and the approximate interfacial geometry (with the cylindrical interfacial surfaces 
replaced by flat planes). 

Figure 14 contains contour plots of the stresses in the epoxy matrix due to a remote imit 
average tensile stress normal to the fiber principal axis. The stress distributions in the glass flbers 
are somewhat featureless. The stresses in the inhomogeneous elements were treated the as the 
homogeneous element stresses in preparing the contour plots. For comparison, the Mmfi 
distribution of matrix stresses is also given in Figure 15 from Reference 5. The latter stresses 
were established using a conventional finite element analysis in which all the elements were 
homogeneous and isotropic. The stress distributions are essentially the same except for a slightly 
higher replacement element stress concentration at the flber/matrix interface along a line of closest 
iqiproach of adjacent fibers in the loading direction. This shows that the replacement scheme can 
give accurate stresses when the stresses and strains within the constituent materials are nonuniform. 
Furthermore, it is not necessary to resort to more refined grids in order to obtain comparable stress 
predictions. 

The transverse Young's modulus prediction from the replacement element solution was 1 .8 
million psi. This also compares favorably with other published values for the same square-packed 
array of glass fibers. For example. Reference 4, lists a value of 1 .7 million psi for a 50% fiber 
volume fraction glass/epoxy with similar constituent properties using conventional finite element 
analyses. 


Sample Problem #4 

This sample problem also represents a 2-D gcnerali 2 ed plane-strain analysis in which the 
constituent material stresses are not uniform. However, the geometry of the reinforcement phase 
was chosen to resemble that of a wavy tow. This microgeometry has sometimes been chosen as 
representative of woven fiber unit cell microgeometries (References 4,6). Figure 16 shows the 
idealized composite structure and a unit cell of that structure. The reinforcing phase consists of 
stacked layers of corrugated aluminum sheets separated by similar layers of epoxy. Perfect 
bonding is assumed between the two phases. The dimensions of the microstructure are given in 
Figure 17. The Young's modulus ofthe composite normal to the plane of Figure 16 can be 
predicted adequately by the rule of mixtures for elements in parallel, but the Young's moduli in the 
X or y- directions require a finite element analysis. This analysis will also consider the 
deformations and stresses in the unit cell as a result of some average strain in the x-direction, with 


13 



all other average strain components held to zero. The constituent material properties ait given in 
Table 1 . The volume fraction of the aluminum is 56%. From symmetry of the microstructuie and 
loading only half of the unit cell needs to be analyzed. 

In order to have a basis of comparison for the approximate analyses a detailed finite 
element analysis was performed on this microstructuie using the NASTRAN code (Ref. 7) and the 
two grids shown in Figure 18. The coarse grid contans 20 elements. The refined grid has 676 
elements. All the elements were homogeneous isotropic CHEXA2 or CWEDGE elements. Three 
independent unit strain cases were run in order to obtain average composite extensional properties 
and the corre^nding stresses and deformations. The average strain case {6 x - 1.0, 6 y ** 6 z 
= 7"yz “ xz * yxy * 0.0} gaw the r^uired jntemal deformations and stresses. The strain 
cases { 6 2 * 1.0, € x= € y*'yyz»yxz*'/xy“0.0} and { € x= € y«* € z* 1,0 
^ y yz= y y xy>°0.0} gave sufficient information to establish the extenslonal 
moduli. The last strain case was obtained by specifying that all average strains vanish and that 
both constituent materials have a unit coefficient of thermal expansion while the unit cell is subject 
to a one degree change in temperature. This was necessary to avoid the occurrence of constant 
displacement terms in the multi-point constraint equations at nodes that were located on surfaces of 
the unit cell where symmetry conditions did not apply (Ref. 7). 

The generalized plane strain, extensional, elastic constants from the NASTRAN models are 



Coarse Grid 

Fine Grid 

Expsi (xlO^) 

3.55 

3.09 

Ey psi (xlO^) 

1.48 

1.43 

Ezpsi (xlO^) 

5.81 

5.83 

Ptx 

0.25 

0.26 

^zy 

0.32 

0.32 

^yx 

0.20 

0.24 


The results from the fine grid are used as the basis of comparison for this example problem. 
Figures 19 and 20 contain plots of the unit cell surface normal deformations and internal stress 
components for the € x ^ 0 strain case. Many of the stress details of the fine grid are not evident 
in the coarse 20-element solution. Even with the refined grid it is not certain whether some of the 
peak stresses have been accurately quantified. The large amount of periodic local bending and 
shearing deformations in the reinforcing sheets are evident in the deformation plots. Large local 
bending stress gradients through the aluminum sheets are also evident in the stress plots. In brief. 
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the task of char^terizing the response of this microstnicture is a more complex problem than the 
iwevious example problem and represents a stiff test of the replacement element method. 

First consider the inhomogeneous element modeling of this microstnicture using the 4 x 4 
grid of Figure 21 and the general energy formulation for the inhomogeneous element stiffness 
matrices. Using four-node, isoparametric, generalized plane strain elements, with the 16-element 
grid the extensional moduli estimates are 

Ex - 4.90 X 10^ psi 

Ey * 3.12 X 10^ psi 

Ex • 5.84 X 10^ psi 

Except for Ez these estimates deviate significantly from the NASTRAN results. If the grid is 
refined from 4x4to8x8as shown in Figure 2 1 the moduli values improve somewhat to 

Ex = 3.95 X 10^ psi 

Ey = 2.32 X 10<5 psi 

Ez * 5.88 X 10^ psi 

However, both the Ex and Ey moduli estimates remain beyond the desired bounds of engineering 
accuracy, and no internal stress data accompany these stiffness estimates. Both of these 
shortcomings can be remedied by the use of replacement elements. 

From the NASTRAN stress results it is obvious that the 4 x 4 grid will not give sufficient 
detail to present any kind of comprehensive picture of the true stress distributions, no matter how 
accurate the replacement clement results may be. Thus the 10 x 10 grid of Figure 21 is applied to 
the current problem with the same type of rectangular replacement element that was used in the 
previous sample problem. With this grid 1 8% of the elements are inhom<^eneous. The resulting 
moduli estimates are 

Ex » 3.21 X 106 psi 

Ey » 1.62 X 106 psi 

Ex = 5.89 X 10^ psi 

These values compare favorably to the base line NASTRAN results. Figure 22 presents the stress 
contours and unit cell surface normal deflections from the 10 x 10 replacement element analysis of 
the € X ^ 0 strain case. The approximations are remarkably consistent with, though slightly less 
detailed than , the fine grid NASTRAN results in Figure 20. The approximations are a major 
improvement in detail over the coarse NASTRAN stress results. 


15 



IV. Three Dimensional Analysis 

The previous sections and example problems have hopefully established the credibility of 
the replacement element method at the one and two-dimensional analysis levels. This section 
extends the method to the 3-D level. Figure 23 shows a parallelepiped element containing two 
different constituent materials. The geometric configuration can be described by specifying the 
volume fraction of one (or both) constituent and the direction of a normal to the interfacial plane. 
The sequence in which the constituent materials appear, as an observer moves along the normal to 
the interfacial plane, is irrelevant to the replacement element method. Figure 23 illustrates the 
two spherical angles \|f , \^2 specify the direction of the nomud to the interfacial plane. These 
two direction angles also serve to locate a set of local coordinates (x,y,z) parallel and normal to the 
interfacial plane. The y and z axes lie in the plane, x is normal to it. The replacement elonent 
concept rearranges the two bulk constituents into a series of parallel plates with the plate surfaces 
paralleling the original interfacial plane. Normal and tangential shear stress continuity is preserved 
across the interface. Compatability of normal strain in the y and z-directions and shear strain in 
the yz plane (of Figiue 23) is maintained across the interfaces. 

Constituent material properties are treated more generally than in the 2-D case. Each 
constituent is assumed to be orthotropic with a plane of isotropy normal to the principal 
reinforcing direction. The principal reinforcing direction must be specified, by means of two 
spherical angles, ((> \ and ((>2. These angles are referenced and measured in the same sense as the 
\|f ] and \|f 2 angles of Figure 23 with the interfacial normal direction replaced by the grain (or 
fiber) direction of the constituent material. Usually the principal reinforcing direction will parallel 
the interfacial plane but this is not assumed in the analysis. 

To form the stress/strain law for the replacement element a number of stress and strain 
transformations must be carried out. Each constituent material has its stress/strain relations 
initially specified in the natural coordinates of the material. These properties must be transformed 
into the x,y,z global coordinates first and then transformed into the x,y,z interfacial coordinates. 
The replacement analysis then yields the replacement material stress/strain law in the x,y,z 
coordinates. Finally, these properties are transformed back into the global x,y,z coordinates for use 
in constructing the element stiffness matrix. This sequence of transformations is retraced (after the 
finite element analysis of the imit cell yields node point deflections and avoage element strains in 
the global coordinates) in order to get constituent material stresses in the natural coordinates of the 
materials. Appendix C derives the replacement element stress/strain equations in the interfacial 
coordinates. Appendix D gives the transformation equations. 
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The 3'D stress and strain transformations are accomplished by a pair of essentially 2-D 
transformations. Each transformation accounts for each ^herical angle of rotation that speciEes 
either the direction of the normal to the interfacial plane or the fiber direction. 


Sample Problem #5 

As in the 1-D and 2-D case, the Hrst 3-D sample problem is an elementry one for which a 
solution is available. However, in this case the blown solution is not exact. The problem concerns 
the "3-D weave" or "XYZ" composite construction (sec Figure 24) in which there are three 
orthogonal Hber directions (Ref 8) . The Ebers remain essentially straight The volume fraction of 
fibers in each of the orthogonal directions usually vary to match the design requirements. The 
types of fibers may also vary with direction. Figure 24 shows one unit cell of the composite 
miciostructure. Symmetry considerations reduce the essential part of the unit cell that must be 
analyzed to one eighth of the total unit cell volume. This reduced volrune is shown in Figure 25. It 
has a 25% volume fraction of interstitial bulk matrix, a 25% volume fraction of unidirectional 
composite with fibers in the x-direction, a 37.5% volume fraction of composite in the y-direction 
and a 12.5% voliune fraction of composite in the z-dircction. The unidirectional material is taken 
to be graphitc/cpoxy with the properties listed in Table 2 under material A. The bulk epoxy 
properties are the same as in the prior sample problems. Using conventional, homogeneous, eight - 
node, isoparametric brick elements and the finite element grid of Figure 26A, the extensional 
composite elastic constants are 


Ex = 5.49 X lO^psi 
Ey = 7.55 X lO^psi 
= 3.43 X lO^psi 
= 0.128 
vi: = O.B. 

\y « 0.055 


The average normal stress in the x-direction in each element as a result of an applied average tensile 
stress of 1000 psi in the global x-direction is given in Figure 27. The results are approximate 
because the stresses are not constant within each brick element. 

The problem can also be addressed using the replacement element approach. For 
example, if the finite element grid of Figure 26B were applied to the XYZ microgeometry there 
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would be three of the eight, equal *sized, brick elements that were inhomogeneous. Using the 
replacement element analysis of Appendix D the pairs of inhomogeneous material in each of these 
three elements can be resolved into three different replacement materials. Using one of these 
replacement materials in each of the inhomogeneous brick elements the finite element analysis can 
proceed as a homogeneous element analysis and the composite stiffnesses and average element 
strains obtained. The same replacement material model may then be used to obtain average 
constituent stresses and strains within each element. These stress predictions are given in Figure 
28. A comparison of Figures 27 and 28 shows that the approximate results from the replacement 
element analysis are of considerable engineering value. The moduli predictions from the two 
models compare as follows: 



Homogeneous Elements 

Replacement Elements 


5.49 X 10^ psi 

5.46 X 10^ psi 

Ey 

7.55 X 10^ psi 

7.55 X 10^ psi 

E. 

3.43 X 10^ psi 

3.44 X 10^ psi 

Vyz 

0.128 

0.128 

Vxz 

0.131 

0.130 

Vxy 

0.055 

0.054 


There are no stiffness discrepancies of any note between the models. The details of the input data 
are given in Appendix E where this sample problem is used to demonstrate the input data 
sequences for the interactive use of the replacement element computer code. 


Sample Problem #6 

The next 3-D sample problem represents a composite comprised of solid glass spheres in an 
epoxy matrix. The volume fraction of the glass reinforcing phase is 25%. The spheres are all the 
same size and are assumed to be packed in a cubic array as shown in Figure 29. The ratio of 
sphere diameter to the spacing distance between centers of adjacent spheres (in the direction of 
closest approach) is 0.684. The problem is the prediction of both the principal Young' modulus in 
the x-direction of Figure 29 and the peak normal matrix stress along the line of closest tqjproach 
of adjacent spheres when the composite has an average remotely applied tensile loading of one psi 
in the x-direction, with all other average stress components equal to zero. 

The problem has no known exact solution but a numerical solution could be obtained with 
any general purpose, 3-D, finite elements code based on the use of conventional, homogeneous. 
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isotropic elements. However, it is of current interest to obtain a solution using rectangular grids 
and replacement elements. 

From symmetry considerations only one octant of a unit cell of the composite needs to be 
analyzed. Figure 30 divides this octant into cubic elements with the 4x4x4 subdivision shown. 
Each element is designated by an iJ4c combination of integers. The i int^er indicates the element 
number along x-axis starting at the origin of Figure 30. j and k are the corresponding element 
counts along the y and z-axes respectively. The 1,1,1 element has one comer on the origin and the 
4,4,4 element is the farthest one from the origin. Table 3 contains the spherical angles (\|f j, \)f2) 
that desig nate the direction of the outward pointing normals from the surface of the glass sphere in 
each element. The table also contains the element volume fractions Ouit are glass and epoxy. This 
is all the input data that is necessary to compute the principal moduli of the composite and the 
stresses in each material of each element using replacement elements. In this example there are 16 
inhomogeneous elements out of a total of 64. Each element is modeled as an 8-node, 
isoparametric, cubic element. The constituent propertes are given in Table 1 . The predicted 
Young's modulus in any of the global coordinate directions of Figure 30 is 0.86x10^ psi. The 
corresponding Poisson's ratio is 0.29 and the shear modlus is 026 x 10^ psi. The peak normal 
stress concentration in the matrix is 2.5. It occurs at the glass/epoxy interface. The stress 
concentration at the same point in a continuous fiber reinforced composite with the same ratio of 
fiber diameter to adjacent fiber spacing is 1 .80. The stresses within the constituent materials of the 
replacement elements appeared to be consistent with the stresses in the neighboring isotropic 
elements. The distribution of normal stress along two faces of the imit cell is shown in Figure 3 1 . 

Sample Problem #7 

The last example of the use of the replacement element analysis considers the plain weave 
unit cell and microgeometry of Figure 32 subjected to uniaxial tension in a reinforcing direction. 

In this model the resin-impregnated and cured tows are considered to be non-circular tubes of 
homogeneous orthotropic material that are woven together. These undulating tubes are bonded 
together at all areas of contact and bonded to the bulk matrix pockets which fill all the interstitial 
gaps between the tubes. The dimensions of the resin filled tows, the tow spacings and the other 
geometric details were chosen to best match the microgeometries observed in photomicrographs of 
woven graphite/epoxy composites (Ref. 1). The analysis was done for the purposes of (a) 
predicting the cxtensional stiffness properties of a thick laminate made from synunetrically stacked 
layers of plain-weave reinforced composite and (b) predicting the detailed stresses and strains 
within one unit cell of this laminate when it is subject to a simple uniaxial tensile stress in one of 
the principal tow reinforcing directions. 
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By the use of structural and load symmetry the essential portion of the unit cell that needs 
to be analyzed can be reduced in volume and complexity. Figure 32 shows one unit cell of the 
plain weave microstructure with four planes of symmetry of both load and structure. Only the 
fraction of unit cell volume between the four planes of symmetry needs to be considered. This 
enclosed volume is shown in Figure 33 with a set of coordinates that parallel the edges of this 
regular hexahedron of essential structure. The origin of the coordinates is at the centroid of the 
hexahedron. These three coordinate axes arc also axes of 180® rotational symmetry of both load 
and structure. Hence, only one quarter of this volume is essential to the analysis. Figure 34 shows 
this reduced volume which represents only one sixteenth of the original unit cell volume. Further 
symmetry exists for the structure but not the loading. Figure 34 also shows a simple retangular 
finite element mesh superposed on the essential structure The use of replacement elements permits 
the application of this grid without much regard for the internal boundaries between the two tow 
materials and the bulk matrix. The mesh has been graded to give added stress detail near the 
crossover point of the upper and lower tows (at the origin of Figure 34). The number of finite 
elements in the smallest essential volume is 64 with 125 node points and 375 degrees of freedom 
prior to the enforcement of the boundary conditions. Examination of the microstreture within each 
finite element shows that six, or 9.4% of these elements, contain all three constituent materials. 

(The two tows are considered to be made from two different materials for bookkeeping 
convenience.) Fourteen, or 2 1 .9% of the elements, contain only one constituent material. The 
remaining 44, or 68.7%, contain two constituent materials. This high percentage of replacement 
elements (78.1%) makes this sample problem different from the previous ones which only required 
a small number of replacement elements. Another es.sential difference is the presrace of elements 
containing three constituents. These special elements are treated as follows. 

First, note that the two tow materials are in direct contact with each other in each element, 
rather than being separated by a layer of bulk matrix. Thus, the reinforced portion of each element 
that contains tow material can be treated as a subelement that contains only two constituent 
materials. Application of the replacement element logic can then be used to combine these two tow 
m at<»riak into a single anisotropic replacement material. One new factor in this reasoning is that 
the subelemcnt containing the two constituents is, in general, no longer a right hexahedron. This 
docs not appear to invalidate the replacement process. After both tow materials have been lumped 
together into a new replacement material then the process can be repeated, combining the new tow 
replacement material with the bulk matrix material. The only new factor in the latter application 
of the method is that the combined tow material may be generally anisotropic. This possibility is 
covered in Appendix C. With these generalizations in place there does not ^pear to be any reason 
to prevent the repeated application of the replacement material logic as many times as necessary in 
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any given clement as long as a "tree diagram" of constituent material combinations, as shown in 
Figure 35 A, can be described. Each outer branch © , @ , (D of the tree diagram represents a 
constituent material. Each junction of two materials © , ® represents an application of the 
replacement material logic. The lower trunk of the tree diagram ® represents the final 
replacement material that is used to form the stiffness matrix for the element. The present analysis 
code (Appendix F) is only general enough to handle the tree diagram of Figure 35A. No more 
complexity was required for this sample problem. 

As examples, consider the microgeometry of a few of the elements from the current sample 
problem. The clement designated 0 in Figure 34 contains only one constituent material, the bulk 
matrix. The tree diagram for this element is a single trunk of one material with no branches or 
junctions. No replacement element analysis is required. 

The clement designated (b| in Figure 34 contains two constituent materials, the bulk 
matrix and one tow material. Figure 35B isolates this clement and shows its tree diagram. The two 
constituent material branches combine at the single junction to form the trunk material. A single 
application of the replacement logic suffices for this element. Figure 35C isolates element 
from Figure 34. This element contains all three constituent materials: the bulk resin and both tow 
materials. Its tree structure is identical to Figure 35A. The replacement logic is applied to the two 
tow materials © and © at junction © initially to form the new material ©. Material ® 
and bulk matrix material © are then combined at junction ® to form the trunk material © 
via the second application of the replacement logic. 

Some comments on the complex mixed boundary conditions on the six surfaces of the plain 
weave structural model are appropriate. Node points on surfaces normal to the z-axis of Figure 36 
have the customary symmetry conditions of zero normal displacements* and zero shear forces. 

The same conditions also apply on the two sides that are at once normal to the x-axis or y-axis but 
not containing cither axis. However, on the two sides containing the coordinate origin the 
rotational symmetry conditions prevail. Node points along either the x or y-axes cannot displace 
normal to the axis and must have a zero applied force component along the axis. A node point 
along either of these two sides (but not on the x or y-axes) must have a corresponding node point 
that is its mirror image on the opposite side of the coordinate axis that is contained within the side 
in which the original node point is located (see Figure 36). The tangential displacements at these 
two imng ft nodcs must be the mirror image of each other (across the intervening coordinate axis). 
The normal displacement must be equal but opposite. The nodal force components normal to the 
side must be mirror images of each other. The nodal force components parallel to the side must be 


* for rigid body and constant strain displacements 
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equal but oppositely directed from their mirror image across the coordinate axis. Along the edges 
of essential structure a combination of the conditions from the intersecting surfaces apply with the 
displacement conditions prevailing over any contradicting force conditions in any specific 
coordinate direction. Along the z-axis of Figure 36 the displacements normal to the z-axis vanish, 
along with the force component parallel to the z*axis. At the coordinate origin all displacements 
vanish. G>mer displacements are determined by the particular strain case being studied, except 
for displacements conditions at comers A,B,C,D of Figure 36. There the aforementioned mixed 
rotational synunetry conditions apply to forces and displacements normal to the faces containing 
thexory-axis. 

Table 4 contains all the geometric information required for each element. These values 
were all obtained by viewing composite photomicrographs and making many sketches of planar 
cuts through the essential structure. It is a chore that would lend itself well to preprocessing. 
However, it is a matter of only a few days work as opposed to the weeks of work associated with 
setting up and checking out a finite element mesh based upon homogeneous elements. 

The tow composite properties used are typical of imidirectional, intermediate modulus, 
graphite/epoxy prepreg. Most prepregs cure out to about 65% fiber volume fraction. The fiber 
volume fraction within a tow of a fabric reinforced composite is generally in the 70% to 75% 
range. This could justify using higher tow composite moduli in the analysis. However, the loss in 
properties due to the weaving process have never been established. The use of the lower properties 
(associated with 65% fiber volume fraction) is an attempt to compensate for fiber breakage, 
misalignment, and other weaving and processing damage. The overall fiber volume fraction for the 
analysis model was 64% with 15% interstitial bulk matrix volume fraction and 85% tow volume 
fraction. The constituent material properties correspond to the epoxy properties of Table 1 and the 
graphitc/epoxy A properties of Table 2. The predicted extensional elastic constants are, with 
reference to the coordinates of Figure 33, 

Ex = Ey = 7.88 X lO^psi 

Ez = 1.69xl06psi 

Vxz “ Vyz = 0.321 

Vxy ** Vyx * 0.048 

As a reference point, the moduli from test data reported in Ref. 1 are 

Ex = 9.13 X lO^psi (warp) 

Ey » 8.83 xlO^psi (fill) 

Vxy = 0.11 
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With conventional laminate theory, for a cross-ply laminate with a 15% thick layer of bulk resin, 
the result would be 

Ex - Ey » 9.68x lO^psi 
Vxy= 0.050 

The stress results are more interesting than the moduli predictions. For a unit remotely 
applied stress in the warp direction, with all other average stresses held at zero, the peak warp tow 
stress has a value of about 4, giving a stress concentration factor of the same amount. This stress 
occurs inside rather than on the surface of the unit cell and away from the cross over point of the 
adjacent tows. It occurs as a result of high bending plus axial strain in the tow that roughly 
parallels the load direction. Figure 37 contains contour plots of the stress in the fiber direction on 
the primary load carrying tow surface. The axial stress in the fill tows are insignificant. Figure 
37 also contains a plot of the axial fiber strain concentration factors based on the ratio of fiber 
longitudinal strain divided by average composite strain in the load direction. These values differ 
significantly from results reported in Ref. 9. The peak fiber strain concentration from the current 
analysis is about 1.5 compared to 2.6 reported in Ref. 9. Also, the location of the peak strains do 
not coincide. The peak strain occurs on the curved portion of the tow surface away from the edges 
of the tow and away from the inflection point of the principal axis of the tow. In Ref. 9 it occtirs at 
tlie edges of the tow at the adjacent tow cross-over point. Plots of the other stress and strain 
components also differ significantly. The two sets of analyses should not be duplicates of each 
other because there were various differences in the models, the constituent properties, the degree of 
mesh refinement, the order of the elements, etc. However, the differences in the results seem larger 
than expected. Differences in tow cross-sectional variation along the tow axis may account for 
much of the discrepancy. In the current analysis very little tow thickness variation was permitted 
because very little was seen in composite photomicrographs. However, in Ref. 9 significant 
necking of the tow thickness (at the sides of the tow) was built into the analysis model near the tow 
crossover point. Some of the strain concentrations could have been the result of these differences 
in cross-sectional modeling. 

In summary, the stress predictions for the sample problem appear to adequately reflect all 
the major combined bending, stretching and shearing effects that were anticipated in the plane 
weave tension analysis. The causes for some of the local strain differences between this analysis 
and that of Ref. 9 remain to be resolved. 

The rotational symmetry boundary conditions that were used with this sample problem are 
not used frequently and were not included in the computer program listed in Appendix F. They 
were used in this problem simply to avoid the necessity of inverting stiffness matrices larger than 
300 square. The program in Appendix F has the more common conditions of geometric unit cell 
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surface synunetry plus load symmetry (and asymmetry) built into it. The same results eould have 
been obtained using the program in Appendix F with some of the larger array dimensions increased 
four fold, and one quarter of the unit cell volume analyzed rather than one sixteenth of the volume. 
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V. Conclusions 


The three-dimensional elastic analysis of complex composite microstructures is made 
difficult by the constraint imposed by conventional finite element analysis on the correspondence 
of internal material interfaces and element boundaries. The concept of a replacement element is 
introduced for the purpose of relaxing this constraint. The replacement element combines the 
constituent materials within an inhomogeneous element into a single anisotropic material to which 
the established finite element procedures may be applied. This constituent material combination 
depends on simple composite mechanics models for parallel bonded plates. This procedure 
involves a physical rearrangement of the materials within the element and therefore represents an 
idealization or approximation of the true material interactions. It has been shown that the use of 
these replacement elements can incur errors on the order of 20% in the predicted stresses within the 
constituents. However, in the more complex problems in which the replacement elements occur 
less frequently the errors in stiffness and internal stress predictions appear to be within a range that 
is acceptable for some engineering applications; namely, trade-off studies that lead to the ranking 
or selection of specific reinforcement microgeometries to meet specific structural requirements. 

Through the use of several example problems of increasing complexity both the application 
and results of the replacement element method are observed. The application is simpler and easier 
than the conventional finite element method in complicated 3-D problems such as those posed by 
many fabric reinforced composite microgeometries. The results are less accurate and less reliable, 
but still acceptable, in view of the statistical variation in unit cell microgeometries and their 
boundary conditions. A large number of finite elements are still required to model a complex 
microstructurc but beyond that point the mechanical analysis is much easier to automate and 
eventually merge with computerized unit cell microgeometry generators, preprocessors and 
postprocessors. The use of replacement elements still requires some skill in the selection of 
rectangular grids which minimize both the number and complexity of the replacement elements. 

It remains to establish guidelines for the use of replacement elements so as to imnimize the 
approximation errors, and also to improve upon the process itself to make it more sensitive to the 
details of the constituent material distribution within an element. The latter tasks could not be 
imdcrtakcn within the seven man-month scope of this effort. 
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Table 1 : Isotropic Constituents 



E X lO^psi 

V 

G X lO^si 

Aluminum 

10.0 

0.30 

4.0 

E Glass 

10.0 

0.25 

4.0 


0.5 

0.35 

0.18 


Table 2: Orthotropic Constituents 



El X lO^psi 

E 2 ,E 3 X 10 ^ psi 

Vl2»Vl3 

V23 

GR/EPA 

18.0 

1.5 

0.23 

0.35 

GR/EPB 

21.0 

1.7 

0.23 

0.30 



Gi 2 » Gi 3 , X 10 ^ psi 

G 23 xlO^psi 

GR/EPA 

0.7 

0.7 

GR/EPB 

0.7 

0.7 
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Table 3: Microgeometry Data for Sample Problem #6 
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Table 3: Microgeometry Data for Sample Problem #6 

(continued) 


Element No. 

W) 


3 . 1.1 

3 . 1.2 

3 . 1.3 

3 . 1.4 

3,2,1 

3 . 2.2 

3 . 2.3 

3 . 2.4 

3 . 3.1 

3 . 3.2 

3 . 3.3 

3 . 3.4 

3 . 4.1 

3 . 4.2 

3 . 4.3 

3 . 4.4 

4 . 1.1 

4 . 1.2 

4 . 1.3 

4 . 1.4 

4 . 2.1 

4 . 2.2 

4 . 2.3 

4 . 2.4 

4 . 3.1 

4 . 3.2 

4 . 3.3 

4 . 3.4 

4 . 4.1 

4 . 4.2 

4 . 4.3 













Table 4: Microgeometry Data for Sample Problem #7 
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Table 4: Microgeometry Data for Sample Problem #7 (continued) 
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Table 4: Microgeometry Data for Sample Problem #7 (continued) 
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Table 4: Microgeometry Data for Sample Problem #7 (continued) 
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Table 5: Resident Constituent Material Properties 
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Figure 3. Sample problem #2 microgeometry. 



Triangular 
F.E. analysis 



Figure 4. F.E. grids for sample problem #2 In principal coordinates. 
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Triangular 
F.E. grid 


Rectangular 
F.E. grid 


Figure 5. F.E. grids for sample problem #2 In global coordinates. 
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Figure 6. Convergence of rectangular Inhomogeneous F.E. solution 
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Figure 7. 2>D Idealization of inhomogeneous element. 
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Figure 8. Unit cell idealization for sample problem #2. 
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Figure 10. Lattice model of unit cell. 
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Figure 11. Unidirectional giass/epoxy composite. 



Figure 12. F.E. grid for unidirectional composite. 
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Figure 16. Microgeometry for sample problem #4. 
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1.07 q 

Figure 17. Microstructural dimensions for sample problem #4. 
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Figure 18. F.E. grids for NASTRAN analysis of sample problem #4. 
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Rgure 20. Baseline internal stress contours from NASTRAN. 
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Figure 25. 3-D weave microstructure. 



Figure 26. Finite element grids for sample problem #7. 
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Figure 30. Glass sphere inclusion microgeometry and grid 
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Figure 31. Surface normal stress in sample problem #6. 



Figure 32. Plain weave unit cell and microgeometry. 
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Figure 35. Tree diagrams for replacement elements. 



52 











First rotation ((|)i) Second rotation ((^ 2 ) 

about the z-axis about the y-axis 

Figure 39. 3-D transformation geometry. 
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Figure 40. F.E. grid description. 
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Appendix A 

(The 2-D Replacement Element Stress/Strain Relation) 


This Appendix considers the replacement of a contiguous volume, filled with two different 
homogeneous isotropic materials, by a single homogeneous orthotropic material. The original 
isotropic materials are assumed to be separated by a single flat interfacial plane. It is also assumed 
that the stresses within each of the original isotropic constituent materials are constant. The pair of 
dissimilar bulk materials are. in effect, replaced by a series of alternating parallel plates with each 
layer having the same elastic properties as one of the isotropic materials. The interfaces between 
each dissimilar plate are parallel to the original bulk constituent interfacial plane. The thicknesses 
of the alternating layers of each material are in the same proportion as the original volume fractions 
of bulk constituent materials. As the plate thicknesses diminish the layered structure may be 
considered to be a composite material in the macroscopic sense. This composite will have a 
principal axis normal to the interfacial plane and the interfacial plane will be a plane of isofropy of 
the composite. The following analysis relates to the generalized plane strain response of this 
composite when the reference plane of the analysis contains the normal to the interfacial plane. 

Consider the dissimilar materials i and j. both of which are homogeneous and isotropic with 
elastic moduli E], vi. Gj and Ej.vj, Gj respectively. The materials are in the form of thin flat sheets 
bonded together to form a composite as shown in Figure 3. Figure 3 shows two umt cells of the 
composite. The volume fractions of the two materials are vj and vj. Coordinate system 1.23 
(shown in Figure 3) has the axes 1 and 3 parallel to the material interfacial plane. Axis 2 is 
normal to that plane. The stress-strain law for material i, for the generalized plane strain case, is 

given by 



A1 



The corresponding equations for materials j are obtained be substituting j for i in equations 
(Al)and(A2). 

From compatibility of displacements, the composite average strains ( E i.E 2.E3. 
Yi 2 ) arc related to the strains in materials i and j as follows: 


€, = e; = 6,i 

62 = Vj 62 + Vj€ 2 * 

63 = € 3 ’ =€ 3 ) 

?.2= ViT;2+Vj7,i. ^ 


(A3) 


From equilibrium and resolution of forces, the composite stresses ( d i,a2,Cl3. ti 2 ) 
are related to the constituent stresses as follows: 


a, =Via,'+Vja,i 
= az‘ = Ua' 
o'3 = Via3*+Vja3^ 
f,2=T|>2 = T,^ • 


(A4) 


This system of 20 equations may be solved for the average composite stress/strain relation 
as follows. First, solve for the constituent stresses using equation (Al) to get 


K' 


1 

> 

B' 

CD 

ZIa 




> = 

B' 

A'' 

b‘ 

< 




B' 

b' 

A' 


i'»‘J 


(A5) 


where A'=Ej(|-J^|)/(l+J^|)n-2l^|) and B* = E{ i^i/( 1+ |-2^|) • 
Since 0 * 2 ' = CTa* * follows that 


B'€|'+A*€2' + B’€3 = B'e,4 A' 62 ' + B'e3* . 
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Thcn,since £| =€,*=€ i' ““<* €3 = €3 = €3' 

B‘e, +a' 62 ’ + b' £3 = b'6 , + A' € 2 ' + B'e 3 . 

This equation may be rewritten as 

(A‘)4-(A')eJ =(B'-B')e,+(B'-B')€3. (A6) 

. j 

Equations (A3) give another relation between ^2^ .namely, 

(Vi)e^+(Vj) €2’ =€2- 
Solving the two foregoing equations for £j and 



Vj(B’-b‘) a' 
V;(B'- 8 ') A' 


Vj(B'-B')] 1^1 

i €2 i (A7) 


where A = A' Vj + aW; • 

Combining the first three lines from (A4) with (A5) gives 
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«Hi ei =[c',d‘,c‘] 


(AIO) 


(a 1 

< €2 * 

wl>w C*=Vj ( - B')/A •"<* 0‘=AyA. 

Combining (A9) and (AIO) gives 
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By interchanging i and j 
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■ 1 0 0- 
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.001. 
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(All) 


(AI2) 


where C^ = Vj( B'-BM/A and = A'/A • 

Inserting equations (A1 1) and (A12) into equation (A8) gives 


' 0*2 


V 




'A'' B'' B*'" 
B'‘ A‘‘ B'' 

gk gk y^k 



0 0 


o'* C‘ 
0 I 



(A 13) 


This represents the composite extensional stress/strain relation. The extensional 
engineering constants for the composite can be obtained from the inverse of the coefficient 
noatrix of equation ( A t3 ) 
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where 


(e 
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B(C-A) 

4^2 «^2 

(B - CD) 

^ ^ «A 
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B(C-A) 

(A-C ) 
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< CTp 





(B-CD) 

B(C-A) 

(AD- B ) 



A= Vi(A' + B'c‘) + Vj(A4B^C^) 

B= VjB'D'+VjB^D^ 

C= ViB'(C' + n +VjB^(c4l) > 

5= ViA'’D' + VjAjD^ 

Oet= A D-2AB +2B^C -CtD 


The clastic constants are 


E,= £ 3 = Det/(a5-B^) 

£ 2 = Det/(A^-C^) 

^=|SS! = ^ = ^ = B(A-C )/Det 
E| t2 t-2 E3 

hi = (c6-b^)/det , 

E| £3 


\ 


> 
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The composite shear stress/strain relation is obtained as follows. 


Since 


T,i=T,|,T,^ = G> 7,2 
(G ' ) Tie ” (G* ) Ti2 = 0 . 



A5 


— (A 14) 


(AI5) 


(A16) 



From (A3) 


‘Vi)7l2 + <Vj>7l'2 = Tl2- 


Solving the two previous equations for 

G' 




V 


T|2 


where V = Vj + Vj G' . 
Since r ,2 = T ,*2 = T 

T|2 = Vi T|2 + Vj Tii 


then 


= (ViGMr4 + (Vj G hylz 
= (V,G')(Gjai2j.KV,GM( G’y» j 

= ^'(Vi+V()7i2 


= iW)y>z- 


(AI7) 


(AI8) 


This is the composite shear stress/strain relation which could also be obtained from the rule of 
mixtures for stiffnesses in series. Equations (A 1 3) and (A 1 8) can be used to establi^ a 
replacement homogeneous stiffness matrix for any such inhomogeneous element. The finite 
element solution provides the nodal displacements and average strains in each element. Equations 
(A3) ,(A7) and (A 1 7)may then be used to obtain the average strains in each constituent material of 
each element. TTie constituent stress/strain laws give the constituent stresses. Failure theories may 
then be applied to each constituent and the interface if desired. 

Note that all the foregoing equations apply even when one of the two materi als is absent. In 
this case the stress and strain predictions for the missing material can be ignored. 
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Appendix B 

(Generalized Plane Strain Transformation of Elastic Constants) 


The transformation equations for plane stress and plane strain are well known (Ref. 4). 
However, for the generalized plane strain case, the recognized equations are incomplete. There is 
another Young’s modulus, another coefficient of mutual influence, and two additional Poisson’s 
ratios related to the out of plane (i.e out of the plane of transformation) response. The added 
Young’s modulus normal to the plane of transformation is invariant with respect to the 
transformation, but the two additional Poisson’s ratios and the extra coefficient of mutual 
influence are not invariants and their transformation equations are not commonly known. This 
Appendix contains their derivation. 


For the coordinate system shown in Figure 26 the stress and strain transformation 
equations, for generalized plane strain, are given by (Ref. 4), 
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1 \ 0 sc / 6 X \ 

^ €2 1 s 2 c2 0 -SC 6y I 

63 I 0 0 1 0 ^ z 

^ 7 i 2 ^ I-2SC 2SC 0 C 2 -S 2 J ^ 7 xy ' 

where S and C designate sin ({> and cos (|) respectively. The stress/strain laws in the 
principal coordinates of an orthotropic material are (Ref 4) 


(B4) 
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- WEa 
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In a coordinate system obtained by a rotation of (}) degrees about the z (or 3) -axis the more 
general anisotropic equations are apply; namely: 

'Ex| r I Ax -VEy -W^z WGxylk’ 

Ey[ -;/xr/Ex I Ay -VEz WCxyAyL-, 

^2 "^xz/Ex -^yzAy • A 2 T^z.xy/Cxy CTz 

A*y] [^xy,x/Ex ^xy.y/Er ^xy-Az 'AxyK- 

From (B2), (B3) and (B5) 



^ o-scl 

0 sc 

t, C| 

0010 -1“-^ i 0 


0 2 SC 
S* 0 - 2 SC 

0 0 r 0 


2 SC- 2 SCOC?q 0 0 0 yi^SCSCO cf-#J Txy 
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Therefore, from (B6) and (B7) after matrix multipUcation 




0- = S' * ( 

Ey E, I 


♦ 


G|2 


2 yi2 
E, 


) 


S^C^+ 


c! 

E2 


* I- f (i* rriJ"''' 

t^3|C2 + Vj2S2 

»'2,= U2,S^*P3zC^ 

2(f +f +4|!?-^)s2e+£(S^+C^) 

6,, \E| Ez E| G|2/ b|2 

J2«y,»= (|.+m2- ' kc^- (f+^-TJ-js^C 

Ex \E| E| G,2/ \E2 E) G|2( 

l7'^= + (ra^TT 

’7x»,2= 2 ( >^32 - »^3|) sc . 


These are the complete transformation equations for the generalized plane strain case. The 
equations preceded by an asterisk are not considered in the plane stress or plane strain case. 

The reciprocal relations, from the required symmetry of the stress/strain coefficient matrix, 
arc given by 

yxy _ ^ = ViM , VMy = 2*ity 

Ex Ey Gxy Ex Gxy Ey 

Ez Ex Ez Ey Gxy Ez 

The last three of these do not appear in the plane stress (or strain) case. 
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Appendix C 

(The 3-D Replacement Element Stress/Strain Relation) 


This Appendix derives the 3-D replacement element stress/strain relation for a subelement 
containing two elastically dissimilar materials. The two constituent materials are generally 
anisotropic and separated by an interfacial plane that parallels the y^ axes of Figure 23. The 
stresses in any constituent material are assumed to be constant throughout the constituent. From 
considerations of equilibrium and action/reaction across the interface the following relations exist 
between the various constituent and average stress components 

a-x = (J-/ = o-x^ 
ay = Vi cTy'+ Vj 

= Vi o-z'-f- Vjai^ 

Ty? = Vi Tyz'+ VjTyl^ 
fxl - Txi '= Txl^ 

Txy= '^xy ' = • 



From geometric compatibility the following constraint conditions apply to the constituent and 
average strains. 

= Vi 6-x' + Vj€xi 

6 -y = 6 -y‘ = 6 -yi 

e-z = e-zV= 

Tyi = - y^^ 

y%z - Vi 7xz'+ Vj yu^ 

7xy = ViTxy'-^VjTxy^ 

Equations (Cl) may be rewritten as 




Cl 



where the matrices [C‘] and [d] are the coefficient matrices of the constituent material stress/strain 
equations in the x,y^ coordinate system of Figure 23. From equations (Cl), (C2) and (C3) 



C2 


where Cnm designates (C^mn * Cinin) C‘mn ts the nm th clement of [C*]. If P] is the inverse 
(obtained numerically) of the leading coefficient matrix in equation (C 4 ) then 


/ • \ 
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\ / 


(C 5 ) 


where T 5 jnn is the mn th element of [D]. Substituting values of the constituent strains from (C 5 ) 
into (C 3 ) gives the following stress/strain law for the replacement material within the element of 
Figure 23 . 
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(C 6 ) 


All=Vi(C|'| D|2+C|5 1^22 + C|g 032)+ V](C|^I D42+C13 D32+ C|g Dg2) 


f r\ 


A 12 “ Vj (C 21 0|2 + C25 P22 + C26 P32) + Vj (C^i D42 + C^5 D52+ 1^02) 


where 



Ai 4 *Vj (C41 D|2 +C45 P22+C4g D32)+Vj (C41 D42+C45 □52**’ ^46 
A15 * Vj ( C51 Di 2 + C55 P22+ C56 D32) + Vj (C51 D42+ C^5 D52+ C56 Dq2 ) 

Vj(Cg| D,2+Cg5D22+CggD32)'*‘Vj(Cg| D42+Cg5D52+CggDg2) 

A22 = V||c 22+C2| ( D|| C21 + B13C25 +D15C26) 

+ C25 ( D21 C21 + ^23^ 25 + D25C26) 

+ C2g( D31 C21 + D33C25 + D35C2g)| 

+ Vj ^21 ^ D41 C21 + D43C25 + D45C26) 

+ C2^ D5, C21 + D53C25 + D55C26) 

+ C2g( Dg| C21 + Dg3C25 + Dg5C26)| 

^23 Vj |C23+ C21 (D|| C31 + D|^ C35 + D|5 C3g) 
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+ Vj |C 23 + C21 (D41 C3, + D43 C35 + D45 C3g) 

+ C2^5 (D51 C31 + D53C35+ D55 €30) 

+ C^g (DgiCj, + Bg3C35+ Dg5C3g)J 

A 24 = Vi|c 24 +C21 (D|| C41 + D13C45 + D|5 C4g) 

+ D21 C4, + D23 C45 + D25 C4g) 

+ C2g(D3|C4| + D33C45+ D35C4g)| 
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0 36 ( D61C3, + O03O35 + 055030)! 
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03s( O21 C41 +D23C45 + D25C40) 

■•'C36( D3, C41+D33C45 + 035040)! 
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+ ^3^5 ( D5, C41 + D53C45 + O55O40) 

+ 03^ ( 001 C41+D03O45 + O05O40)| 


A35® V, (C31 0,4+035 O24+O30 034)+ Vj (O31 0 i| 4 *^^ 35 ^% 4 '^^ 36 ^^ 
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A44“ V||C44+ C4|(D|| C^l ■** ^13^45 ^15 ^46 ) 
+ 045(021 C 41 + D 23 C 45 + D 25 €40 ) 
^46(^31^41 +D 33 C 45 + D 35 C 40 )} 

'*‘yi{^44'*' ^41 (*^4I ^41 ^45^45 ^45 C 46 ^ 

'*'^45^^51^41 ^53^45 **■ ^55 ^46 ^ 

+ C4g(0g,C4| + Dg3C45 + Dg5 €40 )| 


A45“Vi(C4^ 

A46=Vi(C4‘, 

A55= Vi(Csi 


C45 ^ 24 "^ C4I5 D34)+Vj (€4^, D44+ €4^ D54+C4‘'g Dg4) 

D,g+ 04*5 D 26 + C4g D3g) + Vj (C 41 D4g+C4^5 05g+C4^g Dg^) 

0,4+ 05*5 D 24 + C^g 034 )+ Vj (Cs^, D 44 + C 5*5 D 54 + Cs^ 6 . 4 ) 


A56=Vi(C5* 
A66= Vi(Cg'| 


•^16^ C 55 ^26^ C 56 03g)+ Vj (Cs^ O 46 + C 4 O 56 + C 56 Ogg) 

D|6+ C65D26+Cgg03g)+Vj (Cg^i D4g+Cg^5 Djg+C^g Dgg). 


The matrix [A] is symmetric. 

The replacement stress/strain equation (C6) must be transformed into the global x,y^ coordinate 
system for use in forming the stiffness matrix of the replacement element. When the finite element 
solution to the deformations of the unit cell is obtained the average strains in each element can be 
computed. In a replacement element the average strains can be transformed into the x,y^ 
interfacial coordinate system. Equations (C2) and (C3) will then give the average strains in each 
constituent material. These strains can be transformed into the principal axis of the material and 
the constituent stresses computed. Any yield or initial failure criteria can then be applied to the 
constituent materials or the interface. 
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Appendix D 

(The 3-D Stress/Strain Transformations) 


Any 3-D stress/strain coefficient matrix in the x,y^ coordinate system of Figure 39 A can 
be transformed into the global x,y^ coodinate system in two stages. The first stage consists of a 
rotation of <1> i degrees about the z-coordinate axis of Figure 39A. The positive sense of rotaton 
for (]) 1 is clockwise to an observer at the origin looking in the positive z-direction. From the 
equilibrium of triangular wedge elements whose faces are normal to the coordinate axes it can be 
shown that the stresses in both coordinates are related as follows 



where S j * sin (j) i and C j = cos (J> \ . The reverse or inverse transformation is readily obtained 
by substituting - ({>i in place of ()>i in the previous equations giving 



D1 



By definilion the two coefficient matrices in (Dl) and (D2) must be the inverses of each other. The 
strain transformation equations for the same rotation of (j)] d^ree about the z-coordinate are 
derived from geometric considerations alone and given by 


Cf Sf 0 0 0 -S,C, 1 

$2 Cp 0 0 0 S,C, €y 

0 0 I 0 0 0 ^ ^ 

0 0 0 C, S, 0 

0 0 0 -S| C| 0 y^2 

2S,C,-2S,C, 0 0 0 (C,2-Sf)j 


with the reverse transformation 

C2 $2 0 0 0 S,C| • 

$2 c2 0 0 0 -S,C, 

0 0 1 0 0 0 

0 0 0 C| -S| 0 

0 0 0 S, C, 0 

-2S,C,2S,C, 0 0 0 (C2-S2) 

The second stage of the transformation consists of a rotation of <}>2 d^rees about the y- 
coordinatc as shown in Figure 39B. In this case the positive sense of rotation for ())2 is counter 
clockwise to an observer at the origin looking in the positive y-dircction. The relationships 
between the stresses and strains before and after this latter transformation are 
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where $2 “ sin (|>2 and C2 “ cos (J>2. 
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The stresses and strains in the final x,y, 2 -coordinates are thus related to the original stresses 
and strains in the x,y^ coordinates by the following equations obtained by multiplying the 
appropriate pair of coefficient matrices in equations (D3) through (D8) in the correct order. 




DA 



If.^ \ 
t X 




s\ S,S,Cj C,SjCj s.c,4 


€x 



s! 

cf 

0 
(/T 

1 

o 

o 

o 


€y 




sV 

cl -SftC, -C,S^Cj s,c,s^ 


> 

Wi 


as,c,s,-2s,c,sj 

0 C,C, -S,Cj -(CVsl)Sj 


7y2 




5|‘cl-^’ c,(c*-s|) -zsfjS^Cj 


Txz 

XY ^ 


-as,c,Cj 2Sf,Cj 

0 C,S, -S,% (CVs^lC.^ 


[TxY;. 


pl2) 


Note that the coefficient matrix in (Dl 1) is the transpose of the coefficient matrix in (DIO) and the 
coefficient matrices in p9) and (D12) are also transposes of each other. This fact considerably 
simplifies the use of these equations. 
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Appendix E 

(Instructions for Use of the 3-D Computer Program for Resolving Stresses in a Unit Cell) 


llie program REPLACE is a FORTRAN code that analyzes a 3-D unit cell structure for 
the composite moduli and the internal stresses in each material of each element of the finite element 
model of the unit cell. The only element in the program is the eight-node isoparametric brick 
element. The stiffness matrix for this element is formed by numerical integration over eight 
Gaussian integration points. The boundary conditions of structural symmetry are assumed to 
apply on all faces of the unit cell. 

The input data can be considered to be four different data packages. The first package 
serves to establish an ad hoc file of constituent materials and their elastic properties for use in the 
current run. The second data package establishes the hexagonal finite element grid for the unit cell. 
The third data package contains the constituent material distribution information for each bnck 
clement in the finite clement grid. The last data package contains the average applied stresses to be 
used for the detailed internal stress computations. Each data package will be described in the 
foregoing sequence, starting with the constituent materials file. 

The program operates interactively and is largely self-explanatory through the prompt 
messages. The first block of input data serves to establish an ad hoc array of constituent matenals 
and their clastic constants for use in the current run. These constituents can be selected from the 
eight sets of resident materials whose properties correspond to unidirectional high, medium and 
low modulus graphite/epoxy, unidirectional glass/epoxy, bulk aluminum, bulk epoxy, etc. (see 
Table 5). New sets of material properties can be input either by inserting new DATA statements at 
the beginning of the program or by following a sequence of material input prompts. 

The first piece of input data is a single digit integer (NM) that specifies the number of 
constituent materials to appear in the ad hoc list of materials for that run. Not all of the matenals 
on the list need to be used and the same material may appear more than once. Each material is 
presumed to be orthotropic with a plane of isotropy. Thus, six clastic constants suffice to define 
the linear response of the material. Whatever, additional material constants are associated with the 
yield or failure criteria must also be input. The six clastic constants (in the sequence in which they 
arc input and stored) are the Young's modulus in the principal reinforcing direction, the Young s 
modulus in the plane of isotropy, the longitudinal/transvcrse Poisson's ratio, the Poisson’s ratio in 
the plane of isotropy, the longitudinal/transvcrse shear modulus, and the shear modulus in the plane 

of isotropy. A maximum strain failure criteria is currently in the program. It requires four 

additional constants per material. The four input constants associated with this criteria arc the 
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longitudinal tension and compression strain and the transverse tension and compression strain. 

The initial input number of materials (NM) establishes a do loop for filling the materiel property 
(MP) array that contains the ten input constants for each of the (NM) constituent materials. The 
program requests a constituent material data location number for data insertion in the first row of 
the MP array. If a material number between one and eight is specified then the ten material 
constants from that data statement are inserted in the first row of the MP array for the properties of 
what will be known subsequently as constituent material number one. If a data statement number 
greater than eight is called for then the program makes ten queries for each material property to be 
inserted in the first row of the MP array. This sequence is repeated until each row of the material 
property array is filled. The constituent material input sequence also establishes a numbering 
scheme for recalling the constituent material properties. The first row of MP array to be filled is 
henceforth material number one, the second is material number two, etc. 

Following the material property selection is the description of the geometry of the 
rectangular finite element grid to be used in the analysis of the unit cell structure. The first input 
quantity specifies the side length of the unit cell in the x-direction. The second quantity is an 
integer (NBX) that specifies the number of brick elements along one side of the unit cell in the x- 
direction. If NBX is greater than one then the distance of each node point from the origin in the x- 
direction must be specified. This is done by specifying the x-distance from the origin to the 
farthest interior point of each clement that lies along the x-axis of Figure 40, starting with the 
nearest element to the origin and ending with the farthest one. Each distance is designated as a 
percentage of the x-side length of the unit cell. The last distance in the x-dircction is not specified, 
but is assumed to be 100% of the unit cell side length in that direction. The same set of quantities 
arc then specified for the y-dircction of the unit cell grid and then repeated once more for the 2 - 
direction. This establishes the finite element grid. 

It remains to describe the material distribution within each brick element. This is 
done by means of a triple nested do loop starting with the brick element closest to the origin 
of Figure 40. The material distribution within that clement is described in its entirety 
before the inner do loop indexes to the next clement in the plus z-direction for the same 
material distribution data. This inner do-loop continues to index in the z-direction until the 
last element touching the z-axis is fully described. The middle do-loop then indexes to the 
second brick element (from the origin) along the y-axis. The inner do-Ioop once again 
ranges over each element in this second stack of brick elements, requisitioning materials 
information for each in the same sequence. When each z-stacking of brick elements along 
the x"K) face of the unit cell is described then the outer do-loop indexes to the next brick 
element along the x-axis and the two inner do-loops arc restarted. Now consider the way 
the distribution of material in each brick element is described. 
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The material distribution within a brick element must be reducible to a series of 
constituent material junctions as shown in Figure 35 A. Each branch or trunk represents a 
different material with the two branch materials combined to form the trunk material (as 
described in Appendix C). The main trunk of the tree structure represents the single 
material to be used in that element stiffness matrix calculation. 

The tree structure can,in principal, be as complex as necessary as long as each 
junction contains no more than two branches and one trunk. However, in practice, each 
finite element has its tree structure limited to two junctions. No more complexity was 
required for the example problems considered. With this limitation a finite element can not 
contain more than three different constituent materials, as shown in Figure 35C. Each of 
the three outer branches must contain a single constituent material chosen from one of the 
sets of material properties established in the MP array. The description of the outer branch 
consists of the material designation number, corresponding to the MP row number, and the 
pair of spherical coordinate angles ( (j) i , <J>2) that specify the grain or fiber direction, as 
shown in Figure 39. Each junction must also contain a description of the volume fraction 
of each branch and the pair of spherical coordinate angles (\|T 1 , 1 ^ 2 ) specify the 
direction of the normal to the intcrfacial plane separating the two branch materials (see 
Figure 23). Before inputting any unit cell analysis problem the tree structure of each finite 
element should be sketched and labeled as shown in Figure 35. 

In summary, fifteen data values are needed to describe the most general two- 
junction material arrangement in each element. These numbers are prompted and input in 
the following sequence (with reference to Figure 35 A): 

a) material property (MP) array row designation for branch ( 1 ) 

b) spherical angle <j)| for fiber direction in branch (T) 

c) spherical angle (|)2 for fiber direction in branch (I) 

d) material property (MP) array row designation in branch ( 5 ) 

e) spherical angle ()> \ for fiber direction in branch Q) 

f) spherical angle ({>2 for fiber direction in branch Q) 

g) volume fraction of branch Q) material at junction (S) 

h) spherical angle qt], for interfacial normal at junction @ 

i) spherical angle q /2 for interfacial normal at junction 

j) material property (MP) array row designation for branch (|) 
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k) spherical angle ifor fiber direction in branch (3) 
i) spherical angle <^2 direction in branch Q) 

m) volume fraction of branch (|) material at junction 

n) spherical angle \}r j for interfacial normal at junction (§) 

o) spherical angle \|f2 for interfacial normal at junction (§) 

If there is only one junction in an element (or no junctions) less input information is 
needed. For one junction only the first nine inputs ((a)-(i)) are needed. For no junctions 
only the first three inputs ((a)-(c)) are needed. To trigger the correct set of prompts, the 
first piece of input data for any element is the integer 0 through 2 that ^>ecifies the number 
of junctions in the tree structiu^. Then the appropriate set of prompts will automatically 
follow in the foregoing sequence. All angles are to be specified in degrees and decimal 
fractions of a degree. Volume fractions of materials are specified in decimal fractions form 
(0.0 to 1.0.) This completes the description of the material content of each element. The 
stress output for each element is given in the reverse order of the sequence of junction 
descriptors. The stresses are given in the principal axes of the constituent material. Each 
material will also have its minimum margin of safety computed based on a maximum 
strain criteria (with respect to the principal axes of the material). 

The only remaining input is the specification of the six components of the 3-D, 
applied, far-field stresses for which the internal stresses in each material in each element are 
to be computed. 

The 3-D weave (or XYZ material) serves as a simple example for controlling and 
re^nding to the unit cell analysis program input prompts. The composite consists of 
three sets of unidirectional graphite/epoxy tows interspersed as shown in Figure 24. The 
tows in the x,y and z direction are all of a different size. The x-direction tow Ells 25% of 
the unit cell. The y-direcrion tow fills 37.5% of the unit cell. The z-direction tow fills 
12.5% of the unit cell. The remaining 25% of the volume is bulk epoxy. The finite element 
mesh could easily be adjusted so that each element was homogeneous. However, for 
illustrative purposes the mesh will be set up such that there are three inhomogeneous 
elements (out of a total of eight) each containing two different constituent materials. The 
grid is chosen such that all the finite elements have the same dimensions, and there are two 
elements stacked in each coordinate direction. There are only two materials needed; 
unidirectional graphite/epoxy (the first material among the DATA statements) and bulk 
epoxy (the sixth material in the DATA statements). Thus, the first three input integers 



declare that two materials are needed and that they are constituent materials number one 
and six. Since material one was listed first it becomes material number one for the rest of 
the run. Material six hereafter becomes material number two. A printout of the series of 
program prompts and responses are given at the end of this Appendix along with the 
stiffness and stress output. The input is echoed in double parentheses to distinguish it from 
the prompts. The unit cell dimensions are 2.0 units in the x,y and z-directions. The center 
node point along each edge of the unit cell divides the edge into two equal lengths. The 
exploded sketch of Figure 28 shows the sequence in which the eight elements are 
described. Elements ® @ and 0 are inhomogeneous. 

The first clement contains both constituent materials: unidirectional graphite/epoxy 
and bulk epoxy, in equal volumes. The fiber direction angles are <J) and (j)2’^0® for 
material ® . The volume fraction for material one is 0.5. Any direction angles can be 
specified for material two, the isotropic bulk epoxy. In this case <j>i=0O and (^2=0® were 
specified. The interfacial normal has \J / 1 = \p 2 = 0^ as its direction. This information 
fully describes the material content of the first element. 

Element two contains only one constituent material, but half of the fibers are going 
in the y^direction and half are going in the Z“dircction. This can be represented by a single 
junction with both branches made from constituent material number one. One branch has 
fiber direction angles of (J) i = 0° and (|)2= 90°. The other branch has fiber direction angles 
of (J>i= 90® and (})2 = 0®. The volume fraction of branch 0 material is 0.5 and the 
intcrfacial normal has the direction angles \|f [ = 0® and \p 2 = 0®. 

The third element is homogeneous in material one with no junctions. The fiber 
directions are (J>i=(J) 2 = 0 ®. 

The fourth clement differs from the first only in the fiber direction angles. The rest 
of the elements are homogeneous with no junctions. 

The average applied stress is 1 000 psi tension in the x-direction with the other stress 
components equal to zero. 

The output consists of the composite moduli and the stresses in the principal axes of 
each constituent material in each element. Minimum margins of safety are also given. 
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INPOT NO. COMPOSITE MATERIALS NEEDED, NH 
(( 2 )) 

SELECT A MATERIAL NUMBER FROM ONE TO TEN 
(( 1 )> 

SELECT A MATERIAL NUMBER FROM ONE TO TEN 
(( 6 )) 

MATERIAL PROPERTY DATA ECHO 

18000000.00 1500000.00 0.23 0.30 700000.00 700000.00 

O.lOOOE-01 O.lOOOE-01 O.lOOOE-01 O.lOOOE-01 

500000.00 500000.00 0.35 0.35 180000.00 180000.00 

O.lOOOE-01 O.lOOOE-01 O.lOOOE-01 O.lOOOE-01 

INPUT SIDE LENGTH OF UNIT CELL IN X DIR. 

(( 2 . 0 )) 

INPUT NO. SUBCELLS (X DIR.) . IN. UNIT CELL 
(( 2 )) 

INPUT DIST.(t) ORIGIN TO UNIT CELL MODE 2 
((50.0) ) 

INPUT SIDE LENGTH OF UNIT CELL IN Y DIR. 

(( 2 . 0 )) 

INPUT NO. SUBCELLS (Y DIR.) IN UNIT CELL 
(( 2 )) 

INPUT DIST.(%) ORIGIN TO UNIT CELL NODE 2 
((50.0)) 

INPUT SIDE LENGTH OF UNIT CELL IN Z DIR. 

(( 2 . 0 )) 

INPUT NO. SUBCELLS (Z DIR.) IN UNIT CELL 
(( 2 )) 

INPUT DIST.(%) ORIGIN TO UNIT CELL NODE 2 
((50.0)) 

INPUT NUMBER OF JUNCTIONS AT LOCATION 111 
(( 1 )) 

INPUT MATL. NO. 1 AT 1 1 1 

(( 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 
((90.0) ) 

INPUT MATL. NO. 2 AT 1 1 1 

{( 2 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 1ST MATL. VOLUME FRACTION 
(( 0.5)) 

INPUT 1ST INTERFACIAL NORMAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND INTERFACIAL NORMAL ANGLE 

(( 0 . 0 )) 

INPUT NUMBER OF JUNCTIONS AT LOCATION 112 
(( 1 )) 

INPUT MATL. NO. 1 AT 112 
(( 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 
((90.0)) 

INPUT MATL. NO. 2 AT 112 
(( 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 
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((90.0)) 

INPOT 2ND FIBER SPBERZCAl. ANGLE 

(( 0 . 0 )) 

INPOT 1ST MATL. VOLDME FRACTION 

(( 0.5)) 

INPOT 1ST INTERFACIAL NORMAL ANGLE 

(( 0 . 0 )) 

INPOT 2ND INTERFACIAL NORMAL ANGLE 

(( 0 . 0 )) 

INPOT NUMBER OF JUNCTIONS AT LOCATION 121 

(( 0 )) 

SPECIFY THE CURRENT MATL. ID. NO. 

(( 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT NUMBER OF JUNCTIONS AT LOCATION 122 
({ 1 >) 

INPUT MATL. NO. 1 AT 1 2 2 

( ( 2 ) ) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT MATL. NO. 2 AT 122 

(( 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 1ST MATL. VOLDME FRACTION 

(( 0.5)) 

INPUT 1ST INTERFACIAL NORMAL ANGLE 

(( 0.0)) 

INPUT 2ND INTERFACIAL NORMAL ANGLE 

(( 0 . 0 )) 

INPUT NUMBER OF JUNCTIONS AT LOCATION 211 

(( 0 )) 

SPECIFY THE CURRENT MATL. ID. NO. 

(( 2 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

, INPUT NUMBER OF JUNCTIONS AT LOCATION 212 

(( 0 )) 

SPECIFY THE CURRENT MATL. ID. NO. 

(( 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

((90.0)) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT NUMBER OF JUNCTIONS AT LOCATION 221 

(( 0 )) 

SPECIFY THE CURRENT MATL. ID. NO. 

(( 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 

INPUT 2ND FIBER SPHERICAL ANGLE 

(( 0 . 0 )) 
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INPUT NUMBER OF JUNCTIONS AT LOCATION 222 

(( 0 )) 

SPECIFY THE CURRENT MATL. ID. NO.' 

(< 1 )) 

INPUT 1ST FIBER SPHERICAL ANGLE 
((90.0) ) 

INPUT 2ND FIBER SPHERICAL ANGLE 

( ( 0 . 0 ) ) 


ELASTIC CONSTANTS OF THE COMPOSITE 

EX,EY,E2 - 5459114.00 7546501.50 3438649.50 

MUY2,MUX2,MUXY - 0.1276 0.1300 0.0544 

MU2Y,MU2X,MUYX - 0.0581 0.0819 0.0752 


INPUT APPLIED STRESSES IN X,Y,Z COORDINATES 


INPUT X NORMAL STRESS 

(( 1000 . 0 )) 

INPUT y NORMAL STRESS 

(( 0 . 0 )) 

INPUT 2 NORMAL STRESS 

( ( 0 . 0 ) ) 

INPUT Y2 SHEAR STRESS 

(( 0 . 0 )) 

INPUT X2 SHEAR STRESS 

(( 0 . 0 )) 

INPUT XY SHEAR STRESS 

( ( 0 . 0 )) 


STRESSES IN ELEMENT NO. 
MATERIAL NO. 1 
NORMAL 1,2,3 -377.26 

SHEAR 23, 13, 12 -0.49 

MINIMUM MARGIN OF SAFETY IS 


STRESSES IN ELEMENT NO. 
MATERIAL NO. 2 
NORMAL 1,2,3 179.47 

SHEAR 23, 13, 12 0.12 

MINIMUM MARGIN OF SAFETY IS 


STRESSES IN ELEMENT NO. 
MATERIAL NO. 1 
NORMAL 1,2,3 -348.27 

SHEAR 23,13,12 -5.62 

MINIMUM MARGIN OF SAFETY IS 


STRESSES IN ELEMENT NO. 
MATERIAL NO. 1 
NORMAL 1,2,3 -96.79 

SHEAR 23,13,12 -5.46 

MINIMUM MARGIN OF SAFETY IS 


STRESSES IN ELEMENT NO. 
MATERIAL NO. 1 


111 

32.66 179.47 

0.70 0.48 

0.9976 


111 

86.62 81.30 

-0.70 0.49 

0.9759 


112 

66.55 292.88 

-5.46 0.35 

0.9976 


112 

292.88 50.10 

0.35 -5.62 

0.9990 


12 1 


E8 



NORMAL 1»2,3 3337.00 €2.39 47.98 
SHEAR 23/13>12 -0.30 1.68 1.93 
MZNIMOM MARGIN OF SAFETY IS 0.9816 


STRESSES IN ELEMENT NO. 1 2 2 
MATERIAL NO. 2 

NORMAL 1,2,3 190.64 91.97 86.67 
SHEAR 23,13,12 -0.11 1.15 -1.28 
MINIMUM MARGIN OF SAFETY IS 0.9744 


STRESSES IN ELEMENT NO. 122 
MATERIAL NO. 1 

NORMAL 1,2,3 -135.34 190.64 17.88 
SHEAR 23,13,12 -1.15 -0.42 1.28 
MINIMUM MARGIN OF SAFETY IS 0.9990 


STRESSES IN ELEMENT NO. 2 1 1 
MATERIAL NO . 2 

NORMAL 1,2,3 134.54 61.80 55.90 
SHEAR 23,13,12 -0.30 -0.54 0.21 
MINIMUM MARGIN OF SAFETY IS 0.9813 


STRESSES IN ELEMENT NO. 2 1 2 
MATERIAL NO. 1 

NORMAL 1,2,3 -106.32 282.81 49.86 
SHEAR 23,13,12 -5.07 -0.75 -5.87 
MINIMUM MARGIN OF SAFETY IS 0.9990 


STRESSES IN ELEMENT NO. 2 2 1 
MATERIAL NO. 1 

NORMAL 1,2,3 3307.09 61.25 44.10 
SHEAR 23,13,12 0.81 0.53 1.53 
MINIMUM MARGIN OF SAFETY IS 0.9818 


STRESSES IN ELEMENT NO. 222 
MATERIAL NO . 1 

NORMAL 1,2,3 -101.96 275.56 46.96 
SHEAR 23,13,12 -1.68 1.24 2.88 
MINIMUM MARGIN OF SAFETY IS 0.9990 


E9 




Appendix F 

(FORTRAN Program Listing) 


FI 



PROGRAM REPLACE 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


cccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CC PROGRAM TO COMPUTE THE 3-D INTERNAL STRESSES CC 
CC IN A UNIT CELL OF AN INCLUSION ARRAY CC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SPECIFICATION STATEMENTS 

PARAMETER (MM-10, MMM-300, NNN-400) 

REAL MP(8,10),KS(24,24) 

DIMENSION SK{24,24),UVWS(24) 

DIMENSION KA(7),LC(24) 

DIMENSION TS(6),TSS(6),TS1(6),TS2(6),ST1(6),ST2(6) 
DIMENSION SIG(6) ,STN(6) 

DIMENSION PR (6, 6) , RP (6, 6) , T (6, 6) , BIG (6, 6) 

DIMENSION DD(6,6) ,DD1(6,6),D02(6,6) 

DIMENSION SS (6, 6) , SSI (6, 6) , SS2 (6, 6) 

DIMENSION SSS1(6,6),SSS2{6,6) 

DIMENSION BM(6,24),DB(6,24) 

DIMENSION PROP (MM, 10) , DX (MM) , DY (MM) , DZ (MM) 

DIMENSION FDX(MM+1) ,FDY(MM+1) ,FDZ(MM-*-l) 

DIMENSION FB(NNN, 7),FS(MMM,7),FT(6,7),TF(6,7),FTT(6,6) 
DIMENSION UVW(NNN, 6) , VU(NNN) 

DIMENSION MNl (MM, MM, MM) , MN2 (MM, MM, MM) ,MN3 (MM,MM,MM) 
DIMENSION FVl (MM, MM, MM) , FV2 (MM, MM, MM) ,NJC (MM,MM,MM) 
DIMENSION ANG1A(MM,MM,MM) ,ANG2A(MM,MM,MM),ANG3A(MM,MM,MM) 
D IMENS ION ANGIB (MM, MM, MM) , ANG2B (MM, MM, MM) , ANG3B (MM, MM, MM) 
DIMENSION AGN1A(MM, MM,MM) ,AGN2A(MM,MM,MM) 

DIMENSION AGN1B(MM,MM,MM) , AGN2B (MM, MM,MM) 

REAL KB(NNN,NNN),KM(MMM,NNN),KN(MMM,MMM) 


C BUILT IN MATERIAL PROPERTY DATA 
C 


DATA (KA(I),I-1,7)/10,6,6,1,6,0,0/, 


1 

(MP 

(1, 

I), 

,I-1,10)/18 

.E6, 

1.5E6, . 

23, . 

30, 

.7E6, 

.7E6, 

.01, . 

.01, 

.01, 

.01/, 

2 

(MP 

(2, 

I) , 

,I-1.10)/2. 

1E7, 

1.7E6, . 

23,. 

30, 

.7E6, 

.7E6, 

.01, . 

01, . 

01, 

.01/, 

3 

(MP 

(3, 

X) . 

, r- 

■1,10) / 3. 

0E7, 

1.7E6, . 

23,. 

30, 

.7E6, 

.7E6, 

.01, . 

01, . 

01, 

.01/, 

A 

(MP 

(4, 

I) - 

, I- 

■1,10)/1. 

0E7, 

1.SE6, . 

25, . 

35, 

.7E6, 

.7E6, 

.01, . 

01, . 

01, 

.01/, 

S 

(MP 

(5, 

I) . 

, I- 

■1,10) / I. 

2E7, 

1.5E6, . 

25,. 

35, 

.7E6, 

.7E6, 

.01, . 

01, . 

01, 

.01/, 

6 

(MP 

(6, 

I) - 

rl- 

•1,10)/ .5E6, . 

5E6, .35 

,.35 

, .18E6, . 

18E6, 

.01, . 

01, . 

01, 

.01/, 

7 

(MP 

(7, 

I)( 

, I- 

■1,10)/1. 

E7,l 

.E7, .30 

, .30 

, .4E7, .4E7,.01,.01 

, .01 

,.01/, 

6 

(MP 

(8, 

I) . 

. I- 

■1,10)/1. 

E7,l 

.E7, .25 

, .25 

, .4E7, .4E7,.01, .01 

, .01 

, .01/ 


C 

C INITIALIZE VARIABLES 



ISYM-0 



DO 10 I 

-1,MMM 


DO 10 J 

-1,7 

10 

FSd, J) 

-0.0 


DO 12 I 

-1,MMM 


DO 12 J 

-1,MMM 

12 

KN(I, J) 

-0.0 


DO 15 I 

-1,6 


DO 15 J 

-1,7 


TF(I, J) 

-0.0 

15 

FT(I, J) 

-0.0 


WRITE (6 

,9100) 


READ (5, 

9030) NM 


WRITE (6 

,9899) NM 


C MATERIAL PROPERTY DATA INPUT 
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DO IB 

WRITE (6, 9180) 

REAO(5,9030) M 
KRITE(6,9899) M 
IF<M.GT.8) THEN 
WRITE (6, 9120) 

READ(5,9010) PROP (1,1) 

WRITE (6, 9130) 

READ(5,9010) PROP (I, 2) 

WRITE (6, 9140) 

READ(5,9010) PROP(I,3) 

WRITE (6, 9150) 

READ(5,9010) PROP (I, 4) 

WRITE (6, 9160) 

READ (5, 9010) PROP (I, 5) 

WRITE (6, 9170) 

READ(5,9010) PROP (I, 6) 

WRITE (6, 9175) 

READ (5, 9015) PROP (I, 7) 

WRITE (6, 9177) 

READ (5, 9015) PROP (I, 8) 

WRITE(6, 9178) 

READ (5, 9015) PROP (I, 9) 

WRITE (6, 9179) 

READ(5,9015) PROP(I,10) 

END IF 

IF(M.LE.8) THEN 
DO 17 J-1,10 

17 PROPd, J)-MP(M, J) 

END IF 

18 CONTINUE 
WRITE(6,9560) 

WRITE(6,9190) 

DO 19 I-1,NM 

WRITE(6, 9020) (PROP (I, J) , J-1, 6) 

19 WRITE(6,9025) PROP (I, 7) , PROP (1, 8) , PROP (I, 9) , PROP (1, 10) 

READ MESH GEOMETRY 

WRITE(6, 9560) 

WRITE(6, 9440) 

READ (5, 9000) XL 
WRITE(6,9898) XL 
WRITE(6, 9080) 

READ (5, 9030) NBX 
WRITE (6, 9899) NBX 
FDX(l)-0,0 
FDX(NBX+1) -100,0 
IF(NBX.LE.l) GO TO 25 
DO 22 I-1,NBX-1 
WRITE(6,9460) I+l 
READ(5,9000) FDX(I+1) 

22 WRITE(6,9898) FDX(I+1) 

25 CONTINUE 

WRITE (6, 9560) 

WRITE (6, 9450) 

READ (5, 9000) YL 
WRITE (6, 9898) YL 
WRITE (6, 9090) 

READ (5, 9030) NBY 
WRITE (6, 9899) NBY 
FDY(1)-0.0 
FDY(NBY+1)-100.0 
IF(NBY.LE.l) GO TO 35 
DO 30 I-1,NBY-1 
WRITE(6, 9460) I-t-1 
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READ(5,9000) FDY(H-l) 

30 WRITE (6, 9896) TOY {1*1) 

35 CONTINUE 

WRITE (6, 9560) 

WRITE (6, 9455) 

READ (5, 9000) ZL 
WRITE (6, 9898) ZL 
WRITE (6, 9095) 

REAO(5,9030) NBZ 
WRITE (6, 9899) NBZ 
FDZ(l)-0.0 
FDZ(NBZ+1) -100.0 
IF(NBZ.LE.l) GO TO 135 
DO 130 I-1,NBZ-1 
WRITE (6, 9460) I+l 
READ (5, 9000) FDZ(I+1) 

130 WRITE (6, 9898) FDZ(I-t-l) 

135 CONTINUE 

NP- (NBX+1) * (NBY+1) * (NBZ+1) *3 
DO 230 I-1,NP 
DO 230 J-1,6 
230 UVW(I,J)-0.0 
DO 250 I-1,NNN 
DO 240 J-1,7 
240 FB(I,J)-0.0 
DO 250 J-1,NNN 
250 KB<I,J)-0.0 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C BEGIN OUTER DO LOOP OVER THE NO. OF BRICK ELEMENTS C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

DO 2400 I-1,NBX 
DO 2400 J-1,NBY 
DO 2400 K-1,NBZ 

GET ELEMENT DIMENSIONS 

A-FDX(I+1)-FDX(I) 

A-A*XL/100. 

AA-0.5*A 
DX(l)-0.0 
DX(2)-A 

B-FDY(J+1)-FDY(J) 

B-B»YL/100. 

BB-0.5*B 
DY(l)-0.0 
DY(2)-B 

C-FDZ(K+1)-FDZ(K) 

C-C*ZL/100. 

CC-0.5*C 
DZ(1)-0.0 
DZ(2)-C 
VOL-A*B*C 

DO 260 11-1,24 
DO 260 JJ-1,24 
260 KS(II, JJ)-0.0 

INPUT TYPE OF MATERIAL JUNCTION IN THE ELEMENT (0,1, OR 2) 

WRITE (6, 9560) 

WRITE(6,9060) I,J,K 
READ(5,9030) NJC(I,J,K) 

WRITE(6,9899> NJC(I,J,K) 

JNC-NJC(I, J,K) 
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c INPUT MATERIAL TYPE AND FIBER DIRECTION ANGLES (ONLY BRANCH) 

C 

IF(JNC.LT.l) THEN 
WRITE (6, 9320) 

READ(5,9030) MN1(I,J,K) 

MN-MNld, J»K) 

WRITE(6,9899) MN1(I,J,K) 

WRITE (6, 9480) 

READ (5, 9000) A1 
ANGlAd, J,K)-A1 
WRITE (6, 9898) A1 
WRITE (6, 9490) 

READ (5, 9000) A2 
ANG1B(I,J,K)-A2 
WRITE (6, 9898) A2 
CALL TRANS2 (A1,A2,T) 

C GET STRESS-STRAIN MATRIX (SS) IN MATL. COORD. SYSTEM 
C 

CALL GETSS(MM,MN,PROP,SS) 

C GET STRESS-STRAIN MATRIX (SS) IN GLOBAL COORDINATES 
C 

DO 290 11-1,6 
DO 290 JJ-1,6 
SUM-0.0 
DO 280 KK-1,6 

280 SUM-SUM+SS(II,KK) *T(KK, JJ) 

290 PR(II, JJ)-SUM 
C 

DO 300 11-1,6 
DO 300 JJ-1,6 
SUM-0 . 0 
DO 295 KK-1,6 

295 SUM-SUM+T(KK, II) *PR(KK, JJ) 

300 BIG (II, JJ) -SUM 
END IF 

C INPUT MATERIAL TYPE AND FIBER DIRECTION ANGLES (FIRST BRANCH) 
C 

IF(JNC.GE.l) THEN 
WRITE(6,9700) I,J,K 
READ(5,9030) MN1(I,J,K) 

WRITE(6,9899) MN1(I,J,K) 

MN-MNl (I, J,K) 

WRITE(6,9480) 

READ(5,9000) ANG1A(I,J,K) 

WRITE(6,9898) ANG1A(I,J,K) 

A1-ANG1A(I, J,K) 

WRITE(6,9490) 

READ(5,9000) ANG1B(I,J,K) 

WRITE(6,9898) ANG1B(I,J,K) 

A2-ANG1B(I, J,K) 

C GET STRESS-STRAIN MATRIX (SS) IN MATL. COORD. SYSTEM 
C 

CALL GETSS(MM,MN,PROP,SS) 

Q 

C GET STRESS-STRAIN MATRIX (SS) IN GLOBAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 

DO 400 11-1,6 
DO 400 JJ-1,6 
SUM-0 . 0 
DO 350 KK-1,6 
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350 SOM-SUM+SS(H,KK) *T(KK, JJ) 

400 PR(I1,JJ)-S0M 
DO 450 IX-1,6 
DO 450 JJ-1,6 
SOM-0 . 0 
DO 440 KK-1,6 

440 SUM-SUM+T{KK, II)*PR(KK, JJ) 

450 SSI (II, JJ) -SUM 
C 

C INPUT MATERIAL TYPE AND FIBER DIRECTION ANGLES (SECOND BRANCH) 
C 

WRITE(6,9710) I,J,K 
READ(5,9030) MN2(I,J,K) 

WRITE(6,9899) MN2(I,J,K) 

MN-MN2(I, J,K) 

HRITE(6,9480) 

READ (5, 9000) ANG2A(I,J,K) 

WRITE{6, 9898) ANG2A(I,J,K) 

A1-ANG2A(I, J,K) 

WRITE (6, 9490) 

READ (5, 9000) ANG2B(I,J,K) 

NRITE(6,9898) ANG2B(I,J,K) 

A2-ANG2B(I, J,K) 

C 

C GET STRESS-STRAIN MATRIX (SS) IN MATL. COORD. SYSTEM , 

C 

CALL GETSS(MM,MN,PROP,SS) 

C 

C GET STRESS-STRAIN MATRIX (SS) IN GLOBAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 

DO 500 H-1,6 
DO 500 JJ-1,6 
SUM-0 . 0 
DO 490 KK-1,6 

490 SUM-SUM+SS(II,KK) *T(KK, JJ) 

500 PR (II, JJ) -SUM 
DO 550 11-1,6 
DO 550 JJ-1,6 
SUM-0 . 0 
DO 540 KK-1,6 

540 SUM-SUM+T(KK, II) *PR(KK, JJ) 

550 SS2(II, JJ)-SUM 
C 

C INPUT FIRST BRANCH VOL. FRACTION AND INTERFACE NORMAL ANGLES 
C 

WRITE(6, 9720) 

READ (5, 9000) FV1(I,J,K) 

WRITE(6,9898) FV1(I,J,K) 

Vl-FVl (I, J,K) 

V2-1.0-V1 
WRITE (6, 9485) 

READ (5, 9000) AGN1A(I,J,K) 

WRITE (6, 9898) AGN1A(I,J,K) 

A1-AGN1A(I, J,K) 

WRITE (6, 9495) 

READ(5,9000) AGN1B(I,J,K) 

WRITE(6,9898) AGN1B(I,J,K) 

A2-AGN1B(I, J,K) 

C 

C GET STRESS-STRAIN MATRICES (SS) IN INTERFACIAL COORDINATES 
C 

CALL TRANSl (A1,A2,T) 

DO 600 11-1,6 
DO 600 JJ-1,6 
SUM-0 . 0 
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TUM-0 .0 
DO 590 KK-1,6 

SUM-SUM+SSl (II,KK) *T(KK, JJ) 

590 TUM-TOM+SS2(II,KK)*T(KK,JJ) 

RP(1I, JJ)-SUM 
600 PR(II, JJ)-TUM 
DO 650 11-1,6 
DO 650 JJ-1,6 
SOM-0 . 0 
TOM-0 . 0 
DO 640 KK-1,6 
SUM-SOM+T (KK, II) *RP (KK, JJ) 

640 TUM-TUM+T(KK, II) *PR{KK, JJ) 

SSS1(II,JJ)-SUM 
650 SSS2(II,JJ)-TUM 
C 

C DO REPLACEMENT MATERIAL SUBSTITUTION AT FIRST JUNCTION 
C 

CALL GETDD(SSS1,SSS2,V1,V2,DD1,DD2) 

DO 750 11-1,6 
DO 750 JJ-1,6 
DUM-0 . 0 
DO 740 KK-1,6 

740 DUM-DUM+V1*SSS1 (II,KK)*DD1(KK, JJ) +V2*SSS2 (II, KK) *DD2 (KK, JJ) 
750 DD(II, JJ)-DUM 

Q 

C GET REPLACEMENT STRESS-STRAIN MATRIX (SS) IN GLOBAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 

DO 790 11-1,6 
DO 790 JJ-1,6 
SUM-0.0 
DO 780 KK-1,6 

780 SUM-SUM+DD(II,KK) *T(KK, JJ) 

790 PR(II, JJ)-SUM 
DO 850 11-1,6 
DO 850 JJ-1,6 
SUM-0 . 0 
DO 840 KK-1,6 

840 SUM-SUM+T(KK,II)*PR(KK, JJ) 

850 BIGdI, JJ)-SUM 
END IF 
C 

C INPUT MATERIAL TYPE AND FIBER DIRECTION ANGLES (THIRD BRANCH) 

C 

IF(JNC.EQ.2) THEN 
WRITE(6,9715) I,J,K 
READ(5, 9030) MN3(I,J,K) 

WRITE(6,9899) MN3(I,J,K) 

MN-MN3(I, J,K) 

WRITE(6,9480) 

READ (5, 9000) ANG3A(I,J,K) 

WRITE (6, 9898) ANG3A(I,J,K) 

A1-ANG3A(I, J,K) 

WRITE (6, 9490) 

READ (5, 9000) ANG3B(I,J,K) 

WRITE (6, 9898) ANG3B(I,J,K) 

A2-ANG3B(I, J,K) 

C 

C GET STRESS-STRAIN MATRIX (SS) IN MATL. COORD. SYSTEM 
C 

CALL GETSS(MM,MN,PROP,SS) 

C 

C GET STRESS-STRAIN MATRIX (SS) IN GLOBAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 
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DO 1400 11-1,6 
DO 1400 JJ-1,6 
SUM-0 . 0 

DO 1350 KK-1,6 

1350 SUM-SUM+SS(I1,KK) *T(KK, JJ) 

1400 PR<11,JJ)-SUM 
DO 1450 11-1,6 
DO 1450 JJ-1,6 
SOM-0.0 

DO 1440 KK-1,6 

1440 SUH-SUM-)-T(KK, 11)*PR(KK,JJ) 

1450 SSI (11, JJ) -SOM 
DO 1550 11-1,6 
DO 1550 JJ-1,6 
1550 SS2(11, JJ)-B1G(11, JJ) 

C 

C INPUT LAST BRANCH VOL. FRACTION AND INTERFACE NORMAL ANGLES 
C 

KRITE(6, 9720) 

READ (5, 9000) FV2(1,J, K) 

WRITE (6, 9898) FV2(I,J,K) 

V1-FV2(I, J,K) 

V2-1.0-V1 
WRITE (6, 9480) 

READ(5,9000) AGN2A(I,J,K) 

WRITE(6,9898) AGN2A(1,J,K) 

A1-AGN2A(1, J,K) 

WR1TE(6,9490) 

READ (5, 9000) AGN2B(I,J,K) 

WRITE (6, 9898) AGN2B(I,J,K) 

A2-AGN2B(I, J,K) 

C 

C GET STRESS-STRAIN MATRICES (SS) IN INTERFACIAL COORDINATES 
C 

CALL TRANSl (A1,A2,T) 

DO 1600 11-1,6 
DO 1600 JJ-1,6 
SUM-0.0 
TUM-0 . 0 

DO 1590 KK-1,6 
SUM-SUM+SSl (II,KK) *T(KK, JJ) 

1590 TUM-TUM+SS2 (II,KK) *T(KK, JJ) 

RP(II, JJ)-SUM 
1600 PR(II, JJ)-TUM 
DO 1650 11-1,6 
DO 1650 JJ-1,6 
SUM-0.0 
TUM-0. 0 

DO 1640 KK-1,6 
SUM-SUM+T (KK, II) *RP (KK, JJ) 

1640 TUM-TUM+T(KK, II) *PR(KK, JJ) 

SSSl (II, JJ)-SUM 
1650 SSS2(II, JJ)-TUM 
C 

C DO REPLACEMENT MATERIAL SUBSTITUTION AT SECOND JUNCTION 
C 

CALL GETDD(SSS1,SSS2,V1,V2,DD1,DD2) 

DO 1750 11-1,6 
DO 1750 JJ-1,6 
DUM-0 . 0 

DO 1740 KK-1,6 

1740 DUM-DUM+V1*SSS1 (II, KK) *DD1 (KK, JJ) +V2*SSS2 (II, KK) *DD2 (KK, JJ) 
1750 DD(II, JJ)-DUM 
C 

C GET REPLACEMENT STRESS-STRAIN MATRIX (SS) IN GLOBAL COORDINATES 
C 
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CALL TRANS2 (A1,A2,T) 

DO 1790 11-1,6 
DO 1790 JJ-1,6 
SUM-0.0 

DO 1700 KK-1,6 

1780 SUM-SUM+DD(II,KK) *T(KK, JJ) 
1790 PR(II,JJ)-SUM 
DO 1850 11-1,6 
DO 1050 JJ-1,6 
SUM-0.0 

DO 1840 KK-1,6 

1040 SUM-SUM+T(KK, II)*PR(KK, JJ) 
1850 BIG (II, JJ) -SUM 
END IF 

1900 CONTINUE 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C BEGIN INTEGRATION SCHEME TO GET ELEMENT STIFFNESS MATRIX C 
C INNER DO LOOP OVER EACH INTEGRATION POINT BEGINS HERE C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


DO 2000 11-1,2 
DO 2000 JJ-1,2 
DO 2000 KK-1,2 


DO 1920 111-1,6 
DO 1920 JJJ-1,24 
1920 BM(III, JJJ)-0.0 
C 

C SUBCELL GEOMETRY CALCULATIONS 
C 

X-0.57735*AA 
IF(II.EQ.l) X— X 
Y-0.57735*BB 
IF(JJ.EQ.l) Y— Y 
Z-0.57735*CC 
IF(KK.EQ.l) Z— Z 
C 

C DO GAUSSIAN INTEGRATION SCHEME 
C 

CALL GETB{AA,BB,CC,X,Y, Z,BM) 

DO 1930 111-1,6 
DO 1930 JJJ-1,24 
DB(III, JJJ)-0.0 
DO 1930 KKK-1,6 

1930 DB ( 1 1 1 , JJJ) -DB (III, JJJ) +BIG (I II , KKK) *BM (KKK, JJJ) 

DO 1950 111-1,24 
DO 1950 JJJ-1,24 
DO 1950 KKK-1,6 

1 9 50 KS (1 1 1 , J J J) -KS ( 1 1 1 , J J J) +BM ( KKK, 1 1 1 ) *DB (KKK, JJJ) 

2000 CONTINUE 

DO 2040 111-1,24 
DO 2040 JJJ-1,24 

2040 KS (III, JJJ)-KS (III, JJJ)* VOL/8.0 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C END OF INNER DO LOOP OVER INTEGRATION POINTS C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C PUT SMALL STIF. MATRIX (KS) INTO BIG STIF. MATRIX (KB) 

C 

DO 2050 11-1,24 
DO 2050 JJ-1,24 
2050 SK(II,JJ)-KS(II,JJ) 

LC ( 1) - ( (NBZ+1 ) * (NBY+1 ) * (I-l) + (NBZ+1) * ( J-D + (K-1) ) *3+1 
LC(7)-LC(1)+(NBZ+1)*3 
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LC(l3)-LC(l) + (NBZ+l)*{NB1f+l)*3 
LC(19)-LC(13)+(NBZ+1)*3 
DO 2200 KK-2,6 
LC<KK)-I.C{KK-1)+1 
LC (KK+6) -LC (KK+5) +1 
LC (KK+12) -LC (KK+11) +1 
2200 LC(KK+18)-LC(KK+17)+1 
DO 2300 11-1,24 
III-LC(Il) 

DO 2300 JJ-1,24 
JJJ-LC(JJ) 

2300 KB(III, JJJ)-KB(III, JJJ)+SK(II,JJ) 

2400 CONTINUE 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C END OF OUTER DO LOOP ON NO.>-OF ELEMENTS IN UNIT CELL C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C CALC. DISP. VECTORS FOR 6 HOMOGENEOUS UNIT STRAIN CASES 
C 

DO 2420 I-1,NBX-»^1 
DO 2420 J-1,NBY+1 
DO 2420 K-1,NBZ+1 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (NBZ+1) * (J-l) + (K-1) ) *3 
UVW(L+1,1)-FDX(I) *XL/100.0 
UVW(L+2,2)-FDY(J) *YL/100.0 
UVW(L+3,3)-FDZ(K)*ZL/100.0 
UVW(L+2,4)-FDZ(K) *ZL/100.0 
DVW(L+1,5)-FDZ{K) *ZL/100.0 
2420 UVW(L+1,6)-FDY{J) *YL/100.0 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C BEGIN NORMAL STRAIN ANALYSIS C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

GO TO 2460 
2440 ISYM-1 

DO 2445 I-1,MMM 
DO 2445 J-1,MMM 
2445 KN(I,J)-0.0 

DO 2447 I-1,MMM 
DO 2447 J-1,NNN 
2447 KM(I,J)-0,0 

DO 2450 I-1,NP 
DO 2450 J-1,7 
2450 FB(I,J)-0.0 
2460 CONTINUE 
C 

C USE ZERO FORCE CONDITIONS TO ELIMINATE INNER FORCES 
C 

IN-0 

IF((NBX.GT.l) .AND. (NBY.GT.l) .AND. (NBZ.GT.l)) THEN 
DO 2510 1-2, NBX 
DO 2510 J-2,NBY 
DO 2510 K-2,NBZ 

L- ( (NBZ+1 ) * (NBY+1) * (I-l) + (NBZ+1) * ( J-1) + (K-l) ) *3 

DO 2500 M-1,NP 

KM(IN+1,M)-KB(L+1,M) 

KM(IN+2,M)-KB(L+2,M) 

2500 XM(IN+3,M)-KB(L+3,M) 

2510 lN-IN+3 
END IF 
C 

C USE ZERO FORCE B.C.S ON Z-NORMAL FACES 
C 

1F( (NBX.GT.l) .AND. (NBY.GT.l) ) THEN 
DO 2530 1-2, NBX 
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DO 2530 J-2,NBY 

L- ( (NBZ+1 ) * (NBY+1 ) * (I-l ) + (NBZ+1 ) * ( J-1 ) ) *3 

LL-L+NBZ*3 

DO 2525 M-1,NP 

KM(IN+1,M>-KB(L+1,M) 

KM(1N+2,M)-KB(L+2,M) 

KM(IN+3,M)-KB(LL+1,M) 

2525 KM(IN-«-4,M)-KB(LL-»-2,M) 

2530 IN-IN+4 
END IF 

OSE ZERO FORCE B.C.S ON Y-NORMAL FACES 

IF( (ISYM.EQ.O) .AND. (NBX.GT.l) .AND. (NBZ.GT.l)) THEN 
DO 2555 1-2, NBX 
DO 2555 K-2.NBZ 

L- ( (NBZ+1) * (NBY+1) * (1-1) + (K-1) ) *3 
LL-L+ (NBZ+1 ) * (NBY) *3 
DO 2550 M-1,NP 
KM(IN+1,M)-KB(L+1,M) 

KM(1N+2,M)-KB(L+3,M) 

KM(IN+3,M)-KB(LL+1,M) 

2550 KM(IN+4,M)-KB(LL+3,M) 

2555 IN-IN+4 
END IF 

IF( (ISYM.EQ.l) .AND. (NBX.GT.l) .AND. (NBZ.GT.O)) THEN 

DO 2570 1-2, NBX 

DO 2570 K-1, NBZ+1 

L-( (NBZ+1) * (NBY+1) * (I-1)+(K-1) ) *3 

LL-L+(NBZ+1)*(NBY)*3 

DO 2565 M-1,NP 

KM(IN+1,M)-KB(L+2,M) 

2565 KM(IN+2,M)-KB(LL+2,M) 

2570 IN-IN+2 
END IF 

USE ZERO FORCE B.C.S ON X-NORMAL FACES 

IF( (ISYH.EQ.O) .AND. (NBY.GT.l) .AND. (NBZ.GT.l)) THEN 

DO 2585 J-2,NBY 

DO 2585 K-2,NBZ 

L-( (NBZ+1) * (J-1)+(K-1) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

DO 2580 M-1,NP 

KM(IN+1,M)-KB(L+2,M) 

KM(IN+2,M)-KB(L+3,M) 

KM(IN+3,M)-KB(LL+2,M) 

2580 KM(IN+4,M)-KB(LL+3,M) 

2585 IN-lN+4 
END IF 

IF( (ISYM.EQ.l) .AND. (NBY.GT.l) .AND. (NBZ.GT.O)) THEN 

DO 2600 J-2,NBY 

DO 2600 K-1, NBZ+1 

I- ( (NBZ+1 ) * (J-1 ) + (K-1 ) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

DO 2595 M-1,NP 

KM(IN+1,M)-KB(L+1,M) 

2595 KM(IN+2,M)-KB(LL+1,M) 

2600 IN-IN+2 
END IF 

USE ZERO FORCE B.C. ON Z-PARALLEL EDGES 

IF( (ISYM.EQ.O) .AND. (NBZ.GT.l) ) THEN 
DO 2630 K-2,NBZ 
L-(K-l) *3 
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LL-L+ (NBZ+1) *NBY*3 
LLL-L+ (NBZ+1) • (NBY+1) *NBX*3 
LLLL-LL+ (NBZ+1) * (NBY+1) *NBX*3 
DO 2620 M-1,NP 

2620 KM ( IN+1 , M) -KB (L+3 , M) +KB (LL+3, M) +KB (LLL+3, M) +KB (LLU>+3,M) 
2630 IN-IN+1 
END IF 
C 

C USE ZERO FORCE B.C. ON Y-PARALLEL EDGES 
C 

IF( (ISYM.EQ.O) .AND. (NBY.GT.l)) THEN 
DO 2680 J-2,NBY 
L-(NBZ+1) *(J-1)*3 
LL-L+NBZ*3 

LLL-L+ (NBZ+1) * (NBY+1) *NBX*3 
LLLL-LL+ (NBZ+1) * (NBY+1)*NBX*3 
DO 2670 M-1,NP 

2670 KM ( IN+1 , M) -KB (L+2 , M) +KB (LL+2, M) +KB (LLL+2, M) +KB (LLLL+2 , M) 
2680 IN-lN+1 
END IF 
C 

C USE ZERO FORCE B.C. ON X-PARALLEL EDGES 
C 

IF( (ISYM.EQ.O) .AND. (NBX.GT.D) THEN 

DO 2710 1-2, NBX 

L- (NBZ+1) * (NBY+1) * (I-l) *3 

LL-L+NBZ‘3 

LLL-L+ (NBZ+1 ) *NBY*3 

LI,LL-LLL+NBZ*3 

DO 2700 M-1,NP 

2700 KM ( IN+1 , M) -KB (L+1 , M) +KB (LL+1, M) +KB (LLL+1,M) +KB (LLLL+1,M) 
2710 IN-IN+1 
END IF 
C 

IJ-IN 

C 

IF(IJ.LT.l) GO TO 2745 
9085 FORMAT (IH ,6F12.0) 

DO 2730 I-1,IJ 
DO 2730 J-1,7 
2730 FS(l,J)-0.0 
DO 2740 1-1,6 
DO 2735 J-1,IN 
DO 2735 K-1,NP 

2735 FS(J,7)-FS(J,7)+KM(J,K) *UVW(K,I) 

DO 2740 J-1,IJ 
FS(J, I)-FS(J,7) 

2740 FS(J,7)-0.0 
2745 CONTINUE 
INN-IN 
IN-0 
C 

C RETAIN INTERNAL DISPLACEMENTS 
C 

IF( (NBX.GT.D .AND. (NBY.GT.l) .AND. (NBZ.GT.l)) THEN 
DO 2755 1-2, NBX 
DO 2755 J-2,NBY 
DO 2755 K-2,NBZ 

L- ( (NBZ+1) * (NBY+1) * (I-l )+ (NBZ+1) * (J-1)+(K-1) ) *3 
DO 2750 M-1,INN 
KN(M, IN+1 ) -KM (M, L+1) 

KN (M, IN+2 ) -KM (M, L+2 ) 

2750 KN(M, IN+3)-KM(M,L+3) 

2755 IN-IN+3 
END IF 
C 
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C USE DISPLACEMENT B.C.S ON Z-NORMAL FACES 
C 

IF( (NBX.GT.l) .AND. (NBY.GT.D) THEN 
DO 2770 1-2, NBX 
DO 2770 J-2,NBY 

L- < (NBZ+1 ) * (NBY+1 ) * ( I -1 ) + (NBZ+1 ) * { J-1 ) ) *3 

LL-L+NBZ*3 

DO 2765 M-1,INN 

KN (M, IN+1 > -KM (M, L+1 ) 

KN(M, IN+2)-KM(M,L+2) 

KN(M,IN+3)-KM(M,LL+l) 

2765 KN(M, IN+4)-KM(M,LL+2) 

2770 IN-IN+4 
END IF 
C 

C USE DISPLACEMENT B.C.S ON Y-NORMAL FACES 
C 

IF( {ISYM.EQ.O) .AND. (NBX.GT.l) .AND. (NBZ.GT.l)) 
DO 2779 1-2, NBX 
DO 2779 K-2,NBZ 

L- ( (NBZ+1) * (NBY+1 ) * (I-l) + (K-1 ) ) *3 
LL-L+ (NBZ+1) * (NBY) *3 
DO 2777 M-1,INN 
KN(M, IN+1)-KM(M,L+1) 

KN(M, IN+2)-KM(M,L+3) 

KN(M,IN+3)-KM(M,LL+l) 

2777 KN(M, IN+4)-KM(M,LL+3) 

2779 IN-IN+4 
END IF 

IF( (ISYM.EQ.l) .AND. (NBX.GT.l) .AND. (NBZ.GT.O)) 

DO 2706 1-2, NBX 

DO 2786 K-1, NBZ+1 

L-( (NBZ+1) * (NBY+1) * (I-1) + (K-1) ) *3 

LL-L+(NBZ+1)*(NBY)*3 

DO 2784 M-1,INN 

KN(M, IN+l)-KM(M,L+2) 

2784 KN(M, IN+2)-KM(M,LL+2) 

2786 IN-IN+2 
END IF 
C 

C USE DISPLACEMENT B.C.S ON X-NORMAL FACES 
C 

IF( (ISYM.EQ.O) .AND. (NBY.GT.D .AND. (NBZ.GT.l)) 

DO 2792 J-2,NBY 

DO 2792 K-2,NBZ 

L-( (NBZ + 1) + (J-D + (K-1) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

DO 2790 M-1,INN 

KN(M, IN+l)-KM(M,L+2) 

KN(M, IN+2)-KM(M,L+3) 

KN(M, IN+3)-KM(M,LL+2) 

2790 KN(M, IN+4)-KM(M,LL+3) 

2792 IN-lN+4 
END IF 

IF ( (ISYM.EQ.l) .AND. (NBY.GT.D .AND. (NBZ.GT.O)) 

DO 2796 J-2,NBY 

DO 2796 K-1, NBZ+1 

L- ( (NBZ+1 ) * ( J-1 ) + (K-1 ) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

DO 2795 M-1,INN 

KN (M, lN+1 ) -KM (M, L+1 ) 

2795 KN(M, IN+2)-KM(M,LL+l) 

2796 IN-IN+2 
END IF 

C 

C USE DISPLACEMENT B.C. ON Z-PARALLEL EDGES 


THEN 


THEN 


THEN 


THEN 
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IF(dSYM.EQ.O) .AND. (NBZ.GT.l)) THEN 

DO 2605 K-2,N6Z 

L-<K-1)*3 

LL-L+(NBZ+1)*NBY*3 
LLL-L+{NBZ+1) * (NBY+1) *NBX*3 
LLLL-LL+ (NBZ+1) * (NBY+1) *N8X*3 
DO 2800 M-1,INN 

2600 KN (M, IN+1 ) “KM (M, L+3 ) +KM (M, LL+3) +KM (M, LLL+3) +KM(M, LLLL+3) 
2605 IN-IN+1 
END IF 
C 

C USE DISPLACEMENT B.C. ON Y-PARALLEL EDGES 
C 

IF((ISYM.EQ.O) .AND. (NBY.GT.l)) THEN 
DO 2840 J“2,NBY 
L-(NBZ+1)*(J-1)*3 
LL-L+NBZ+3 

LLL-L+ (NBZ+1) * (NBY+1) *NBX*3 
LLLL-LL+ (NBZ+1) * (NBY+1) *NBX*3 
DO 2830 M-1,INN 

2830 KN (M, IN+1 ) -KM (M, L+2 ) +KM (M, LL+2) +KM (M, LLL+2) +KM (M, LLLL+2 ) 
2840 IN-IN+1 
END IF 
C 

C USE DISPLACEMENT B.C. ON X-PARALLEL EDGES 
C 

IF( (ISYM.EQ.O) .AND. (NBX.GT.l) ) THEN 

DO 2680 1-2, NBX 

L-(NBZ+1) * (NBY+1) * (I-l) *3 

LL-L+NBZ*3 

LLL-L+ (NBZ+1) *NBY*3 

LLLL-LLL+NBZ*3 

DO 2870 M-1,INN 

2870 KN (M, IN+1) -KM (M, L+1 ) +KM (M, LL+1) +KM(M, LLL+1) +KM(M, LLLL+1) 
2880 IN-IN+1 
END IF 
C 

C GET UNCONSTRAINED DISPLACEMENTS FOR UNIT STRAIN CASES 
C 

CALL MATINV(KN,MMM, IJ,FS,7,7,DET) 

C 

C GET A COMPLETE SET OF TOTAL DISPLACEMENTS 
C 

IN-0 

C 

C ON INTERIOR 
C 

IF( (NBX.GT.l) .AND. (NBY.GT.l) .AND. (NBZ.GT.l) ) THEN 
DO 2905 1-2, NBX 
DO 2905 J-2,NBY 
DO 2905 K-2,NBZ 

L-( (NBZ+1) * (NBY+1) * ( I-l )+ (NBZ+1) * (J-1)+(K-1) ) *3 
IF(ISYM.EQ.l) GO TO 2902 
DO 2900 M-1,3 

UVW(L+1,M)-UVW(L+1,M) -FS(IN+1,M) 

UVW (L+2 , M) -UVW ( L+2 , M) -FS ( IN+2 , M) 

2900 UVW(L+3,M)-UVW(L+3,M)-FS(IN+3,M) 

GO TO 2905 
2902 CONTINUE 
M-6 

UVW (L+1, M) -UVW (L+1, M)-FS (IN+1, M) 

UVW (L+2, M) -UVW (L+2, M) -FS(IN+2,M) 

UVW (L+3 , M) “UVW (L+3 , M) -FS ( IN+3 , M) 

2905 IN-IN+3 
END IF 
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c 

C ON Z-NORMAL FACES 
C 

IF( (NBX.GT.l) .AND. (NBY.GT.D) THEN 
DO 2930 1-2, NBX 
DO 2930 J-2,NBY 

L-( (NBZ+1) * (NBY+1) * (I-1) + (NBZ+D* (J-1) ) *3 
LL-L+NBZ*3 

IF(ISYM.EQ.l) GO TO 2925 
DO 2920 M-1,3 

0VW(L+1,M)-UVW(L+1,M)-FS(IN+1,M) 

UVW (L+2 , M) -UVW ( L+2 , M) -FS ( IN+2 , M) 

OVW (LL+1,M) -UVW (LL+1, M) -FS (IN+3,M) 

2920 OVW (LL+2,M) -OVW (LL+2, M> -FS (IN+4,M) 

GO TO 2930 
2925 CONTINUE 
M-6 

UVW(L+1,M)-UVW(L+1,M) -FS(IN+1,M) 
UVW(L+2,M)-UVW(L+2,M) -FS(IN+2,M) 

OVW (LL+1, M) -UVW (LL+1, M) -FS(IN+3,M) 

UVW (LL+2, M) -UVW (LL+2, M) -FS (IN+4 ,M) 

2930 IN-IN+4 
END IF 
C 

C ON Y -NORMAL FACES 
C 

,IF( (ISYM.EQ.O) .AND. (NBX.GT.l) .AND. (NBZ.GT.l) ) THEN 
DO 2950 1-2, NBX 
DO 2950 K-2,NBZ 

L-( (NBZ+1) * (NBY+1) * (I-l) + (K-1) ) *3 
LL-L+ (NBZ+1) * (NBY) *3 
DO 2945 M-1,3 

UVW (L+1 , M) -UVW ( L+ 1 , M) -FS ( IN+1 , M) 
UVW(L+3,M)-OVW(L+3,M) -FS(IN+2,M) 

UVW (LL+1, M) -UVW (LL+1, M) -FS(IN+3,M) 

2945 UVW(LL+3,M)-OVW(LL+3,M) -FS(IN+4,M) 

2950 IN-IN+4 
END IF 

IF( (ISYM.EQ.l) .AND. (NBX.GT.l) .AND. (NBZ.GT.O)) THEN 

DO 2960 1-2, NBX 

DO 2960 K-1, NBZ+1 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (K-1) ) *3 

LL-L+(NBZ+1)*(NBY)*3 

M-6 

UVW (L+2, M) -UVW (L+2, M)-FS( IN+1, M) 

UVW (LL+2, M) -OVW (LL+2, M)-FS (IN+2, M) 

2960 IN-IN+2 
END IF 
C 

C ON X-NORMAL FACES 
C 

IF( (ISYM.EQ.O) .AND. (NBY.GT.l) .AND. (NBZ.GT.l)) THEM 

DO 2980 J-2,NBY 

DO 2980 K-2,NBZ 

L-( (NBZ+1) * (J-1) + (K-1) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

DO 2975 M-1,3 

UVW (L+2, M) -UVW (L+2, M)-FS (IN+1, M) 

OVW (L+3,M) -UVW (L+3,M)-FS( IN+2, M) 
UVW(LL+2,M)-UVW(LL+2,M) -FS(IN+3,M) 

2975 UVW (LL+3, M) -UVW (LL+3,M) -FS (IN+4,M) 

2980 IN-IN+4 
END IF 

IF( (ISYM.EQ.l) .AND. (NBY.GT.l) .AND. (NBZ.GT.O)) THEM 
DO 2996 J-2,NBY 
DO 2996 K-1, NBZ+1 
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L-((NBZ+1)*(J-1)+(K-1))*3 

( (HBZ+1) * (NBY+1) *HBX) *3 

M*6 

0VW(L+1,M)-0VW(L+1,M)-FS(IH+1#M) 

UVW (LL+1, M) -OVW (LL+1, M) -FS (IH+2,M) 

2996 IN-IN+2 
END IF 
C 

C ON Z-PARALLEL EDGES 

^ IF( (ISYM.EQ.O) .AND. (NBZ.GT.l) ) THEN 

DO 3010 K-2,NBZ 
L-(K-l) *3 

LL-L+(NBZ+1)*NBY*3 
LLL-L+(NBZ+1)*(NBY+1) *NBX*3 
LLLL-LL+ (NBZ+1) * (NBY+1) *NBX*3 
DO 3000 M-1,3 

UVW(L+3,M)-UVW(L+3,M) -FS(IN+1,M) 

UVW (LL'f3, H) -OVW (LL+3,M) -FS (IN+1,M) 
OVW(LLL+3,M)-UVW(LLL+3,M) -FS(IN+1,M) 

3000 UVW (LLLL+3, M) -UVW (LLLL+3, M) -FS (IN+1,M) 

3010 lN-IN+1 
END IF 
C 

c ON Y-PABALLEL EDGES 
C 

IF( (ISYM.EQ.O) .AND. (NBY.GT.l)) THEN 
DO 3040 J-2,NBY 
L-(NBZ+1)*(J-1)*3 
LL-L+NBZ*3 

LLL-L+ (NBZ+1) * (NBY+1) *NBX*3 
LLLL-LL+ (NBZ+1) * (NBY+1) *NBX*3 
DO 3030 M-1,3 

UVW(L+2,M)-UVW(L+2,M) -FS(IN+1,M) 

UVW (LL+2, M) -UVW (LL+2, M) -FS (IN+1,M) 

OVW ( LLL+2 , M ) -UVW ( LLL+2 , M) -FS ( IN+1 , M) 

3030 OVW (LLLL+2,M) -UVW (LLLL+2,M) -FS (IN+1,M) 

3040 IN-IN+1 
END IF 
C 

c ON X-PARALLEL EDGES 
C 

IF( (ISYM.EQ.O) .AND. (NBX.GT.l) ) THEN 

DO 3070 1-2, NBX 

L- (NBZ+1) * (NBY+1) *(I-1) *3 

LL-L+NBZ*3 

1.LL-L+ (NBZ+1) *NBY*3 

LLLL-LLL+NBZ‘3 

DO 3060 M-1,3 

UVW (L+1,M)-UVW(L+1,M)-FS (IN+1, M) 

UVW (LL+1 , M) -UVW (LL+1, M) -FS (IN+1,M) 
UVW(LLL+1,M)-UVW(LLL+1,M) -FS(IN+1,M) 

3060 UVW ( LLLL+ 1 , M) -UVW ( LLLL+1 , M) -FS ( IN+1 , M) 

3070 IN-IN+1 
END IF 

IF (ISYM.EQ.O) GO TO 2440 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccc 

C BEGIN ((GAMMA-YZ) - 1.0) SHEAR STRAIN ANALYSIS C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC(XC 

c 

DO 3145 I-1,MMM 
DO 3145 J-1,MMM 
3145 KN(I,J)-0.0 

DO 3147 I-1,MMM 
DO 3147 J-1,NNN 
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3147 KM(l/J)-0.0 

DO 3150 1-1, NP 
DO 3150 J-1,7 
3150 FB(I,J)-0.0 

C OSE ZERO FORCE CONDITIONS TO ELIMINATE INNER FORCES 
C 

IF ( (NBX.LE. 1) .OR. (NBY.LE. 1) .OR. (NBZ .LE.l) ) GO TO 3165 
DO 3160 1-2, NBX 
DO 3160 J-2,NBY 

DO 3160 K-2,NBZ . . . 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (NB2+1) • (J-1) + (K-1) ) *3 

DO 3155 M-1,NP 
KM(IN+1,M)-KB(L+1,M) 

KM(IN+2,M)-KB(L+2,M) 

3155 KM(IN+3,M)-KB(L+3,M) 

3160 IN-IN+3 
3165 CONTINUE 
C 

c OSE ZERO FORCE B.C.S ON Z -NORMAL FACES 

^ IF( (NBX.LE. 0) .OR. (NBY.LE. 1)) (»3 TO 3240 

DO 3230 I-1,NBX+1 
DO 3230 J-2,NBY 

L-( (NBZ+D* (NBY+1) * (I-l) + (NBZ+1) * (J-1) ) *3 

LL-L+NBZ*3 

DO 3220 M-1,NP 

KM ( IN+1 , M) -KB (L+3 , M) 

3220 KM(IN+2,M)-KB(LL+3,M) 

3230 IN-IN+2 
3240 CONTINUE 
C 

c USE ZERO FORCE B.C.S ON Y-NORMAL FACES 
C 

IF( (NBX.LE. 0) .OR. (NBZ. LE.l) )GO TO 3260 
DO 3255 I-1,NBX+1 
DO 3255 K-2,NBZ 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (K-1) ) *3 
LL-L+(NBZ+1)*(NBY)*3 
DO 3250 M-1,NP 
KM(IN+1,H)-KB(L+2,M) 

3250 KM(IN+2,M)-KB(LL+2,M) 

3255 IN-IN+2 
3260 CONTINUE 
C 

c USE ZERO FORCE B.C.S ON X-NORMAL FACES 
C 

IF( (NBY.LE. 1) .OR. (NBZ. LE.l)) GO TO 3290 

DO 3280 J-2,NBY 

DO 3280 K-2,NBZ 

L-((NBZ+1)*(J-1)+(K-1))*3 

LL-L+ ( (NBZ+l) * (NBY+1) *NBX) *3 

DO 3270 M-l,NP 

KM(IN+1,M)-KB(L+2,M) 

KM ( IN+2 , M) -KB (L+3 , M) 

KM ( IN+3 , M) -KB ( LL+2 , M) 

3270 KM(IN+4,M)-KB(LL+3,M) 

3280 IN-IN+4 
3290 CONTINUE 
3320 CONTINUE 
C 

IJ-IN 

C 

C GET UNCONSTRAINED NODAL FORCES 
C 
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ir<IJ.LT.l) GO TO 3345 

DO 3330 

J-7 

3330 rS(I,J)-0.0 

1-4 

DO 3335 J-1,IN 
DO 3335 K-1,NP 

3335 rS ( J, 7) -FS ( J, 7) +KM( J, K) *UVW (K, I) 

DO 3340 

rS(J, I)-FS(J,7) 

3340 rS(J,7)-0.0 
3345 CONTINUE 
C 

C RETAIN INTERNAL DISPLACEMENTS 
C 

INN-IN 

IN-0 

IF(IJ.LT.l) GO TO 3595 

IF( (NBX.LE.l) .OR. (NBT.LE.l) .OR. (NBZ.LE.l) ) GOTO 3400 
DO 3390 1-2, NBX 
DO 3390 J-2,NBY 
DO 3390 K-2,NBZ 

L- ( (NBZ+1 ) * {NBY+1 ) * (I-l ) + (NBZ+1) * ( J-1) + (K-1) ) *3 
DO 3380 M-1,INN 
KN(M, IN+1)-KM(M,L+1> 

KN(M, IN+2)-KM(M,L+2) 

3380 ICN(M, IN+3)-KM(M,L-*-3) 

3390 lN-IN+3 
3400 CONTINUE 
C 

C USE DISPLACEMENT B.C.S ON Z-NORMAL FACES 
C 

IF( (NBX.LE.O) .OR. (NBY.LE.l) ) GOTO 3440 
DO 3430 I-l,NBX->-l 
DO 3430 J-2,NBY 

L-( (NBZ+1) * (NBY+1) * (I-D + (NBZ+1) * (J-1) ) *3 

LL-L+NBZ+3 

DO 3420 M-1,INN 

KN(M, IN+l)-KM(M,L+3) 

3420 KN(M, IN+2)-KM(M,LL+3) 

3430 IN-IN+2 
3440 CONTINUE 
C 

C USE DISPLACEMENT B.C.S ON Y-NORMAL FACES 
C 

IF( (NBX.LE.O) .OR. (NB2.LE.1) ) GO TO 3480 
DO 3475 I-1,NBX+1 
DO 3475 K-2,NBZ 

L- ( (NBZ+1 ) * (NBY+1) * (I-l) + (K-1) ) *3 
LL-L+ (NBZ+1) * (NBY) *3 
DO 3470 M-1,INN 
KN(M, IN+l)-KM(M,L+2) 

3470 KN(M, IN+2)-KM(M,LL+2) 

3475 IN-IN+2 
3480 CONTINUE 
C 

C USE DISPLACEMENT B.C.S ON X-NORMAL FACES 
C 

IF( (NBY.LE.l) .OR. (NBZ.LE.D) GO TO 3494 

DO 3492 J-2,NBY 

DO 3492 K-2,NBZ 

L-( (NBZ+1) * (J-1) + (K-1) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

DO 3490 M-1,INN 

KN (M, IN+1 ) -KM (M, L+2 ) 

KN (M, IN+2 ) -KM (M, L+3 ) 
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KN(M,IN+3>-KM(M,LL+2) 

3490 KN(M, IN+4)-KM(M,LL+3) 

3492 IN-IN+4 
3494 CONTINUE 
C 

3595 CONTINUE 

GET UNCONSTRAINED DISPLACEMENTS FOR UNIT STRAIN CASES 


CALL MATINV(KN,MMM, IJ,FS,7,7,DET) 

GET A COMPLETE SET OF TOTAL DISPLACEMENTS 
ON INTERIOR 

IF( (NBX.LE.l) .OR. (NBT .LE.l) .OR. (NBZ.LE.l) ) GO TO 3610 

DO 3605 1-2, NBX 
DO 3605 J-2,NBY 

DO 3605 K-2,NBZ . 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (NBZ+1) * (J-1) + {K"1) ) 3 

M"*4 

UVW(L+1,M)-UVW(L+1,M) -FS(IN+1,M) 

UVW(L+2,M)-UVW(L+2,M) -FS(IN+2,M) 

UVW(L+3,M)-UVW(L+3,M) -FS(IN+3,M) 

3605 IN-IN+3 
3610 CONTINUE 
C 

C ON Z-NORMAL FACES 

Q 

IF((NBX.LE.0) .OR. (NBY.LE.D) GOTO 3640 
DO 3630 I-1,NBX+1 
DO 3630 J-2,NBY 

L- ( (NBZ+1) * (NBY+1) * (I-D + (NBZ+1) * (J-D ) *3 
LL-L+NBZ*3 

UVW(L+3,M)-UVW(L+3,M) -FS (IN+1,M) 

UVW (LL+3,M) -UVW (LL+3,M) -FS (IN+2,M) 

3630 IN-IN+2 
3640 CONTINUE 
C 

C ON Y-NORMAL FACES 

Q 

IF ( (NBX.LE.O) .OR. (NBZ.LE. 1) ) <50 TO 3690 
DO 3680 I-1,NBX+1 
DO 3680 K-2,NBZ 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (K-1) ) *3 

LL-L+ (NBZ+1 ) * (NBY ) * 3 

M-4 

UVW(L+2,M)-UVW(L+2,M)-FS(IN+1,M) 

UVW (LL+2 , M) -UVW (LL+2 , M) -FS ( IN+2 ,M) 

3680 IN-IN+2 
3690 CONTINUE 
C 

c ON X-NORMAL FACES 

Q 

IF( (NBY.LE.D .OR. (NBZ.LE. D) GOTO 3790 

DO 3780 J-2,NBY 

DO 3780 K-2,NBZ 

L- ( (NBZ+1 ) * ( J-1 ) + (K-1 ) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

M-4 

UVW ( L+2 , M) -UVW ( L+2 , M) -FS ( IN+1 , M) 

UVW ( L+ 3 , M ) -UVW ( L+ 3 , M) -FS ( IN+2 , M) 

UVW (LL+2, M) -UVW (LL+2, M) -FS (IN+3,M) 
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3775 UVW (LL-«-3, M) -UVW (LL+3, M) -FS (IN+4 ,M) 

3780 IN-IN+4 
3790 CONTINUE 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccc 
C BEGIN ((GAMMA-XZ) - 1.0) SHEAR STRAIN ANALYSIS C 
cccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

DO 3845 I-1,MWM 
DO 3845 

3845 KN(I,J)-0.0 

DO 3847 I-1,MMM 
DO 3847 J-1,NNN 
3847 KM{I,J)-0.0 

DO 3850 I-1,NP 
DO 3850 J-1,7 
3850 FB(I,J)-0.0 
C 

C USE ZERO FORCE CONDITIONS TO ELIMINATE INNER FORCES 
C 

IN-0 

IF( (NBX.LE.l) .OR. (NBY.LE.l) .OR. (NBZ.LE.l) ) GOTO 3865 
DO 3860 1-2, NBX 
DO 3860 J-2,NBY 
DO 3860 K-2,NBZ 

L- ( (NBZ+1 ) * (NBY+1 ) * <1-1 ) + (NBZ+1) * ( J-1 ) + (K-1) ) *3 

DO 3855 M-1,NP 

KM(IN+1,M)-KB(L+1,M) 

KM(IN+2,M)-KB(L+2,M) 

3855 KM(IN+3,M)-KB(L+3,M) 

3860 IN-IN+3 
3865 CONTINUE 
C 

C USE ZERO FORCE B.C.S ON Z-NORMAL FACES 
C 

IF( (NBX.LE.l) .OR. (NBY.LE.O)) GO TO 3940 
DO 3930 1-2, NBX 
DO 3930 J-1, NBY+1 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (NBZ+1) + (J-1) ) *3 

LL-L+NBZ*3 

DO 3920 M-1,NP 

KM(IN+1,M)-KB(L+3,M) 

3920 KM(IN+2,M)-KB(LL+3,M) 

3930 IN-IN+2 
3940 CONTINUE 
C 

C USE ZERO FORCE B.C.S ON Y -NORMAL FACES 
C 

IF( (NBX.LE.l) .OR. (NBZ.LE.l) ) GO TO 4050 
DO 3980 1-2, NBX 
DO 3980 K-2,NBZ 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (K-1) ) *3 
LL-L+(NBZ+D* (NBY)*3 
DO 3970 M-1,NP 
KM(IN+1,M)-KB(L+1,M) 

XM(1N+2,M)-KB(L+3,M) 

KM(IN+3,M)-KB(LL+1,M) 

3970 KM(IN+4,M)-KB(LL+3,M) 

3980 IN-IN+4 
4050 CONTINUE 
C 

C USE ZERO FORCE B.C.S ON X-NORMAL FACES 
C 

IF( (NBY.LE.O) .OR. (NBZ.LE.l) ) GOTO 4150 
DO 4080 J-1, NBY+1 
DO 4080 K-2,NBZ 
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I- ( (NB2+1 ) * ( J-1 ) + (K-1 ) ) *3 
LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 
DO 4070 

KM(IN+1,M)-KB{L+1,M) 

4070 KM(IN+2,M)-KB(LL+1,M) 

4080 IN>lN-»^2 
4150 CONTINUE 
4320 CONTINUE 
C 

IJ-IN 

C 

C GET UNCONSTRAINED NODAL FORCES 
C 

IF(IJ.LT.l) GO TO 4345 

DO 4330 I-1,IJ 

J-7 

4330 FS(I,J)-0.0 
1-5 

DO 4335 J-1,IN 
DO 4335 K-1,NP 

4335 FS ( J, 7) -FS ( J, 7) +KM ( J, K) *UVW (K, I» 

DO 4340 J-1,IJ 
FS(J,I)-FS(J,7) 

4340 FS(J,7)-0.0 
4345 CONTINUE 
C 

C RETAIN INTERNAL DISPLACEMENTS 
C 

INN-IN 

IN-0 

IF(IJ.LT.l) GO TO 4595 

IF( (NBX.LE.l) .OR. (NBY.LE.l) .OR. (NBZ.LE.D) GO TO 4400 
DO 4390 1-2, NBX 
DO 4390 J-2,NBY 
DO 4390 K-2,NBZ 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (NBZ+1) • (J-1) + (K-1) > *3 
DO 4380 M-1,INN 
KN(M, IN+1)-KM(M,L+1) 

KN(M, IN+2)-KM(M,L+2) 

4380 KN(M, IN+3)-KM(M,L+3) 

4390 IN-IN+3 
4400 CONTINUE 
C 

C USE DISPLACEMENT B.C.S ON Z-NORMAL FACES 
C 

IF( (NBX.LE.l) .OR. (NBY.LE.O) ) GO TO 4440 
DO 4430 1-2, NBX 
DO 4430 J-1,NBY+1 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (NBZ+1) * (J-1) ) *3 

LL-L+NBZ*3 

DO 4420 M-1,INN 

KN (M, IN+1 ) -KM (M, L+3 ) 

4420 KN(M, IN+2)-KM(M,LL+3) 

4430 IN-IN+2 
4440 CONTINUE 
C 

C USE DISPLACEMENT B.C.S ON Y-NORMAL FACES 
C 

IF( (NBX.LE.l) .OR. (NBZ.LE.l) ) GOTO 4480 
DO 4475 1-2, NBX 
DO 4475 K-2,NBZ 

L-( (NBZ+1 )* (NBY+1 )* (I-l) + (K-1) ) *3 
LL-L+ (NBZ+1) * (NBY) *3 
DO 4470 M-1,INN 
KN(M, IN+D-KM (M,L+1) 

KN (M, IN+2 ) -KM (M, L+3 ) 
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KN(M, IN+3>-KM(M,LL+l) 

4470 KN(M,INf4)-KM(M,I.I.+3) 

4475 lN-IN+4 
4400 CONTINUE 
C 

C USE DISPLACEMENT B.C.S ON X-NORMAL FACES 
C 

IF({NBY,LE.O) .OR. (NBZ.LE,!)) GO TO 4494 

DO 4492 J-1,NBY+1 

DO 4492 K-2,NBZ 

L- ( (NBZ+1 ) * ( J-1 ) + {K-1 ) ) *3 

LL-L+ ( (NBZ+1) * (NBY+1) *NBX) *3 

DO 4490 M-1,INN 

KN(M, IN+1>-KM(M,L+1) 

4490 KN(M, IN+2)-KM(M,LL+l) 

4492 IN-IN+2 
4494 CONTINUE 
C 

4595 CONTINUE 
C 

C GET UNCONSTRAINED DISPLACEMENTS FOR UNIT STRAIN CASES 
C 

IJ-IN 

CALL MATINV(KN,MMM, IJ,FS,7,7,DET) 

C 

C GET A COMPLETE SET OF TOTAL DISPLACEMENTS 
C 

C ON INTERIOR 
C 

IN-0 

IF( (NBX.LE.l) .OR. (NBY.LE.l) .OR. (NBZ.LE.l)) GO TO 4610 
DO 4605 1-2, NBX 
DO 4605 J-2,NBY 
DO 4605 K-2,NBZ 

L- ( (NBZ+1 ) * (NBY+1) *(!-!)+ (NBZ+1) * (J-1) + (K-1) ) *3 
M-5 

UVW(L+1,M)-UVW(L+1,M) -FS(IN+1,M) 
UVW(L+2,M)-UVW(L+2,M)-FS(IN+2,M) 
UVW(L+3,M)-UVW(L+3,M)-FS(IN+3,M) 

4605 IN-IN+3 
4610 CONTINUE 
C 

C ON Z-NORMAL FACES 
C 

IF { (NBX.LE.l) .OR. (NBY.LE.O) ) GO TO 4640 
DO 4630 1-2, NBX 
DO 4630 J-1, NBY+1 

L- ( (NBZ+1 ) • (NBY+1 ) * (I-l) + (NBZ+1) * (J-1) ) *3 

LL-L+NBZ*3 

M-5 

UVW(L+3,M)-UVW(L+3,M)-FS(IN+1,M) 

UVW (LL+3,M) -UVW (LL+3,M) -FS (IN+2,M) 

4630 IN-IN+2 
4640 CONTINUE 
C 

C ON Y-NORMAL FACES 
C 

IF( (NBX.LE.l) .OR. (NBZ.LE.l)) GO TO 4690 
DO 4680 1-2, NBX 
DO 4680 K-2,NBZ 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (K-1) ) *3 

LL-L+ (NBZ+1) * (NBY) *3 

M-5 

UVW(L+1,M)-UVM(L+1,M)-FS(IN+1,M) 

UVW(L+3,M)-UVW(L+3,M)-FS(IN+2,M) 

UVW (LL+1 , M) -UVW (LL+1 , M) -FS ( IN+3, M) 
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OVW (LL+3, M) -UVW (LL+3, M) -PS (XN4-4,M) 

4680 IN-IN<»-4 
4690 CONTINUE 
C 

C ON X-NORHAL PACES 
C 

IP((NBY.LE.O) .OR. (NBZ.LE.D) GOTO 4790 

DO 4780 J-1,NBY+1 

DO 4780 K-2,NBZ 

L-{ (NBZ+1) * (J-1) + (K-1) ) *3 

LL-L+ { (NBZ+1) * (NBY+1) *NBX) *3 

M-5 

DVW(L+l,M)-OVW(L+l,M)-PS(IN+l,M) 

UVW (LL+ 1 , M) -UVW (LL+1 , M) -PS ( lN+2 , M) 

4780 IN-IN+2 
4790 CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C COMPUTE NODAL PORCES AND ELASTIC CONSTANTS C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

IF(NP.LT.l) GO TO 4950 
DO 4900 I-1,NP 
DO 4900 J-1, 6 
PB(I, J)-0.0 
DO 4900 K-1,NP 

4900 PB ( I, J) -FB ( I , J) +KB { I , K) *UVW (K, J) 

4950 CONTINUE 
C 

C COMPUTE SIDE LOADS FOR EACH UNIT STRAIN CASE 
C 

XNA-YL*ZL 

YNA-XL*ZL 

ZNA-XL*YL 

DO 5400 M-1,6 

DO 5250 J-1, NBY+1 

DO 5250 K-1, NBZ+1 

L-((NBZ+1)*(J-1)+(K-1))*3 

FT(1,M>-FT(1,M)-FB(L+1,M) 

FT ( 6 , M) -FT ( 6 , M) -FB { L+ 2 , M) 

5250 FT(5,M)-FT(5,M)-FB(L+3,M) 

FT(1,M)-FT(1,M) /XNA 

FT(5,M)-FT(5,M) /XNA 

FT(6,M)-FT(6,M) /XNA 

DO 5300 I-1,NBX+1 

DO 5300 K-1, NBZ+1 

L-( (NBZ+1) * (NBY+1) * (I-D + (K-1) ) *3 

FT(2,M)-FT(2,M)-FB(L+2,M) 

5300 FT(4,M)-FT(4,M)-FB(L+3,M) 

FT(2,M)-FT(2,M) /YNA 
FT(4,M)-FT(4,M)/YNA 
DO 5350 I-1,NBX+1 
DO 5350 J-1, NBY+1 

L- ( (NBZ+1) * (NBY+1) * (I-l) + (NBZ+1) * (J-1) ) *3 
5350 FT(3,M)-FT(3,M)-FB(L+3,M) 

FT(3,M)-FT(3,M)/ZNA 
5400 CONTINUE 

DO 5420 1-1,6 
DO 5420 J-1, 6 
5420 FTT(I,J)-FT(I, J) 

C 

C CALCULATE THE ELASTIC CONSTANTS OF THE UNIT CELL 
C 

CALL INV(FTT) 

DO 5450 1-1,6 
DO 5450 J-1, 6 
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o n o o o o 


5450 FT(I,J)-FTT(1,J) 

EX-1.0/FT(1,1) 

EY-1.0/FT(2,2) 

EZ-1.0/FT(3,3) 

GXY-1.0/FT(6,6) 

GYZ-1.0/FT(4,4) 

GXZ-1.0/FT(5,5) 

PRXY— FT(2,1)/FT(1,1) 

PBYX— FT(1,2)/FT(2,2) 

PRXZ— FT(3,1)/FT(1,1) 

PBZX— FT(1,3)/FT(3,3) 

PRYZ— FT(3,2)/FT{2,2) 

PRZY— FT(2,3)/FT(3,3) 

C MIXXY-FT<4,1)/FT(1,1) 

C MIXXZ-FT(6,l)/FTa,l) 

C MIXYZ-FT(5,1)/FT(1,1) 

C MIYXY-FT(4,2)/FT(2,2) 

C MIYXZ-FT(6,2)/FT(2,2) 

C MIYYZ-FT(5,2)/FT{2,2> 

C MIZXY-FT(4,3)/FT(3,3) 

C MIZXZ-FT(6,3) /FT(3, 3) 

C MIZYZ-FT{5, 3) /FT(3, 3) 

C CXYXZ-FT{4, 5) /FT{5, 5) 

C CXYYZ-FT(4, 6)/FT(6, 6) 

C CXZYZ-FT(5, 6) /FT(6, 6) 

WRITE (6, 9560) 

WRITE (6, 9560) 

WRITE (6, 9600) 

WRITE (6, 9560) 

WRITE (6, 9500) EX,EY,EZ 
WRITE(6,9510) GYZ,GXZ,GXY 
WRITE (6, 9520) PRYZ, PRXZ, PRXY 
WRITE(6,9525) PRZY,PRZX,PRYX 
WRITE (6, 9530) MIXYZ,MIYYZ,MIZYZ 
WRITE (6, 9540) MIXXZ,MIYXZ,MIZXZ 
WRITE (6, 9550) MIXXY,MIYXY,MIZXY 
WRITE (6, 9560) 

WRITE (6, 9560) 

READ THE APPLIED STRESSES 

WRITE (6, 9610) 

WRITE(6,9560) 

WRITE(6,9560) 

WRITE(6,9620) 

READ (5, 9010) SX 
WRITE (6, 9895) SX 
WRITE (6, 9630) 

READ{5,9010) SY 
WRITE(6,9895) SY 
WRITE (6, 9640) 

READ (5, 9010) SZ 
WRITE (6, 9895) SZ 
WRITE(6,9650) 

READ (5, 9010) SYZ 
WRITE (6, 9895) SYZ 
WRITE (6, 9660) 

READ(5,9010) SXZ 
WRITE (6, 9895) SXZ 
WRITE (6, 9670) 

READ(5,9010) SXY 
WRITE (6, 9895) SXY 
WRITE (6, 9560) 

WRITE (6, 9560) 

C CALCULATE THE AVERAGE STRAINS 
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STN(1)-SX*FT{1,1)+SY*FT(1,2)+SZ*FT(1,3)+SY2TT{1,4) 

STN (1 ) -STN { 1 ) +SXZ*FT (1, 5) +SXY*FT (1, 6) 

STN (2) -SX*FT (2, 1) +SY*FT (2, 2) +SZ*FT (2, 3) +SYZ*FT(2, 4) 
STN(2)-STN(2)+SXZ*FT(2,5)+SXY*FT(2,6) 

STN (3) -SX*FT (3, 1) +SY*FT (3, 2) +SZ*FT (3, 3) +SYZ*FT(3, 4) 
STN(3)-STN(3)+SXZ*FT{3,5)+SXY*FT(3, 6) 
STN(4)-SX*FT(4,1)+SY*FT(4,2)+SZ*FT(4,3)+SYZ*FT(4,4) 
STN(4)-STN(4)+SXZ*FT(4,5)+SXY*FT(4,6) 
STN(5)-SX*FT(5,1)+SY*FT(5,2)+SZ*FT(5,3)+SYZ*FT(5,4) 

STN (5) -STN (5) +SXZ*FT (5, 5) +SXY*TT(5, 6) 

STN (6) -SX*FT (6, 1) +SY*FT (6, 2) +SZ*FT(6, 3) +SYZ*FT(6, 4) 
STN(6)-STN(6)+SXZ*FT(6,5)+SXY*FT(6,6) 

C 

C GET NODAL DISPLACEMENTS CORRESPONDING TO THE AVG. STRAINS 
C 

DO 5500 I-1,NP 

VU(I)-STN(1)*UVW(I, l)+STN(2)*OVW(I,2)+STN(3)*OVW(I,3) 

VU (I) -VU (I) +STN (4 ) *UVW(I, 4) +STN (5) *OVW{I, 5) +STN (6) *OVH(I, 6) 
5500 CONTINUE 
C 

cccccccccccccccccccccccccccccccccccccccccccccc 
C BEGIN DO LOOP FOR STRESSES IN EACH ELEMENT C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C ESTABLISH ELEMENT DIMENSIONS 

c 

DO 8000 I-1,NBX 
DO 8000 J-1,NBY 
DO 8000 K-1,NBZ 
A-FDX(I+1)-FDX(I) 

A-A*XL/100. 

AA-0.5»A 

B-FDY(J+1)-FDY(J) 

B-B*YL/100. 

BB-0.5*B 

C-FDZ(K+1)-FDZ(K) 

C-C*ZL/100. 

CC-0.5*C 
DO 5600 111-1,6 
DO 5600 JJJ-1,24 
5600 BMdII, JJJ)-0.0 
X-0.0 
Y-0.0 
Z-0.0 

c 

c GET STRAIN / DISPLACEMENT MATRIX 
C 

CALL GETB(AA,BB,CC,X,Y,Z,BM) 

C 

C GET CORNER DISPLACEMENTS 
C 

L- ( <NBY+1) * (NBZ+1 ) * (I-l) + (NBZ+1) * ( J-1) + (K-1) ) *3 

LL-L+(NBZ+1) *3 

LLL-L+ (NBZ+1) * (NBY+1) *3 

LLLL-LLL+(NBZ+1)*3 

DO 5700 M-1,6 

DVWS(M)-VU(L+M) 

UVWS(M+6)-VU(LL+M) 

UVWS (M+12) -VU (LLL+M) 

5700 UVWS(M+18)-VU(LLLL+M) 

C 

C GET AVERAGE ELEMENT STRAINS 
C 

DO 5800 M-1,6 
TS(M)-0.0 
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DO 5B00 N-1,24 

5800 TS(M)-TS(M)+BM(M,N)*OVWS(N) 

DUM-TS{4) 

TS(4)-TS(5) 

TS(5)-TS(6) 

TS(6)-D0M 

C 

c RECALL NUMBER OF JUNCTIONS IN ELEMENT 
C 

JTC-NJC(I, J,K) 

C RECALL MATL. NO. AND FIBER DIRECTION ANGLES (ONLY BRANCH) 
C 

IF{JTC.LT.l) THEN 
MN-MNld, J,K) 

A1-ANG1A(I, J,K) 

A2-ANG1B(I, J,K) 

C 

C GET AVERAGE STRAINS IN MATERIAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T> 

DO 5900 M-1,6 
SIG(M)-0.0 
DO 5900 N-1,6 
SS(M,N)-0.0 

5900 SIG(M)-SIG(M)+T{M,N)*TS(N) 

C WRITE(6,9031) I,J,K 

C WRITE (6, 9020) (SIG (M) ,M-1, 6) 

C 

C GET STRESSES IN THE MATERIAL 
C 

CALL GETSS(MM,MN,PROP,SS) 

CALL GETMS(MM,MN, PROP, SIG, SAFE) 

DO 6010 M-1,6 
TSl (M)-O.O 
DO 6010 N-1,6 

6010 TSl (M)-TS1(M)+SS(M,N)*SIG(N) 

WRITE (6, 9560) 

WRITE(6,9690) I,J,K 
WRITE (6, 9692) MN 

WRITE (6, 9694) TSl (1) ,TS1 (2) ,TS1 (3) 

WRITE (6, 9696) TSl (4) , TSl (5) , TSl (6) 

WRITE (6, 9698) SAFE 
WRITE (6, 9560) 

END IF 
C 

C RECALL INTERFACIAL NORMAL DIRECTION ANGLES (FIRST BRANCH) 
C 

IF(JTC.EQ.l) THEN 
A1-AGN1A(I, J,K) 

A2-AGN1B(I, J,K) 

C 

C GET AVG. STRAINS IN INTERFACIAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 

DO 6050 11-1,6 
DUM-0 . 0 

DO 6040 JJ-1,6 
6040 DUM-DUM+T(II, JJ)*TS(JJ) 

6050 TSS(II)-DUM 
C 

C RECALL MATL. NO. AND FIBER DIRECTION ANGLES (FIRST BRANCH) 
C 

MN-MNl (I, J,K) 

A1-ANG1A(I, J,K) 

A2-ANG1B(I, J,K) 
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c 

C GET STRESS/ STRAIN MATRIX IN GLOBAL COORDINATES 
C 

CALL GETSS(MM,MN,PROP,SS) 

CALL TRANS2(A1,A2,T) 

DO 6400 11-1,6 
DO 6400 JJ-1,6 
SOH-0 . 0 

DO 6350 KK-1,6 

6350 SDM-SUM-^SS(II,KK)*T(KK,JJ) 

6400 PR(II, JJ)-S0M 
DO 6450 11-1,6 
DO 6450 JJ-1,6 
SUM-0 . 0 

DO 6440 KK-1,6 

6440 SOM-SUM+T(KK, II)*PR(KK,JJ) 

6450 SSKII, JJ)-SUM 
C 

C RECALL MATL.NO. AND FIBER DIRECTION ANGLES (SECOND BRANCH) 

C 

MN-MN2(I, J,K) 

A1-ANG2A(I, J,K) 

A2-ANG2B(I, J, K) 

C 

C GET STRESS/STRAIN MATRIX IN GLOBAL COORDINATES 
C 

CALL GETSS(MM,MN,PROP,SS) 

CALL TRANS2(A1,A2,T) 

DO 6500 11-1,6 
DO 6500 JJ-1,6 
SUM-0.0 

DO 6490 KK-1,6 

6490 SUM-SUM-t-SS(II,KX)*T(KK, JJ) 

6500 PR(II, JJ)-SUM 
DO 6550 11-1,6 
DO 6550 JJ-1,6 
SUM-0 . 0 

DO 6540 KK-1,6 

6540 SUM-SUM-^T(KK, II)*PR(KK, JJ) 

6550 SS2(II, JJ)-SUM 
C 

C RECALL FIRST BRANCH VOL.FRRCT. AND INTERFACIAL NORMAL ANGLES 
C 

V1-FV1{I, J,K) 

V2-1.0-V1 
A1-AGN1A(I, J,K) 

A2-AGN1B(I, J,K) 

C 

C GET S.S. MATRICES IN INTERFACIAL COORDINATES 
C 

CALL TRANSl (A1,A2,T) 

DO 6600 11-1,6 
DO 6600 JJ-1,6 
SUM-0 . 0 
TUM-0 . 0 

DO 6590 KK-1,6 
SUM-SUM-^SSl (II, KK) *T(KK, JJ) 

6590 TUM-TUM-«-SS2(II,KK)*T(KK,JJ) 

RP(II, JJ)-SUM 
6600 PR(II, JJ)-TUM 
DO 6650 11-1,6 
DO 6650 JJ-1,6 
SUM-0 . 0 
TUM-0 . 0 

DO 6640 KK-1,6 
SUM-SUM+T (KK, II ) *RP (KK, JJ) 
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6640 TUM“T0M-*-T (KK, II) *PR(KK» JJ) 

SSSKII. JJ)“SOM 
6650 SSS2(II,JJ)-TUM 

C DO REPLACEMENT MATERIAL ANALYSIS AT FIRST JUNCTION 

^ CALL GETDD(SSS1,SSS2,V1,V2,DD1,DD2) 


C 

C GET 
C 


6711 


CONSTITUENT STRAINS IN INTERFACIAL COORDINATES 


DO 6711 11-1,6 
TSl(II)-0.0 
TS2(II)-0.0 
DO 6711 JJ-1,6 
TSKID-TSKII) 
TS2 (II)-TS2 (II) 


+DD1(II, JJ)*TSS(JJ) 
+DD2(II, JJ)*TSS(JJ) 


C GET CONSTITUENT STRAINS IN GLOBAL COORDINATES 
C 

DO 6720 11-1,6 
DUM-0 . 0 
TUM-0 . 0 

DO 6715 JJ-1,6 
DUM-DOM+T (II, JJ> *TS1 (JJ) 

6715 TUM-TUM+T (II, JJ) *TS2 ( JJ) 

TS(II)-DUM 
6720 TSS(II)-TOM 


C GET CONSTITUENT STRAINS IN MATERIAL COORDINATES 
C 

Al— ANGIA ( I, J, K) 

A2-ANG1B(I, J,K) 

CALL TRANS2 (Al,A2,T) 

DO 6730 11-1,6 
DUM-0 . 0 

DO 6725 JJ-1,6 
6725 DUM-DUM+T(II, JJ)*TS(JJ) 

6730 TSKID-DUM 

A1-ANG2A(I, J,K) 

A2-ANG2B(I, J,K) 

CALL TRANS2 (A1,A2,T) 

DO 6740 11-1,6 
DUM-0 . 0 

DO 6735 JJ-1,6 

6735 DUM-DUM+T(II, JJ)*TSS(JJ) 

6740 TS2(II)-DUM 


C 

C 

C 


GET CONSTITUENT STRESSES IN MATERIAL COORDINATES 

MN-MNKI, J,K) 

CALL GETSS(MM,MN,PROP,SS) 

CALL GETMS (MM, MN, PROP, TSl, SAFE) 

DO 6750 11-1,6 
STl(II)-0.0 
DO 6750 JJ-1,6 

6750 STl (II) -STl (II) +SS (II, JJ) *TS1 (JJ) 

WRITE (6, 9560) 

WRITE(6,9690) I,J,K 
WRITE (6, 9692) MN 

WRITE (6, 9694) STl (1) , STl (2) , STl (3) 

WRITE (6, 9696) STl (4) , STl (5) ,ST1 (6) 

WRITE (6, 9698) SAFE 
WRITE (6, 9560) 

MN-MN2 (I, J,K) 

CALL GETSS(MM,MN,PROP,SS) 

CALL GETMS (MM,MN, PROP, TS2, SAFE) 
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OO 6760 11-1,6 
ST2(II)-0.0 
DO 6760 JJ-1,6 

6760 ST2 (I I ) -ST2 ( I I ) +SS ( II , J J) *TS2 ( JJ) 

WRITE (6, 9560) 

WRITE(6, 9690) I,J,K 
WRITE (6, 9692) MN 

WRITE (6, 9694) ST2 (1) , ST2 (2) ,ST2 (3) 

WRITE (6, 9696) ST2 (4) , ST2 (5) ,ST2 (6) 

WRITE (6, 9698) SAFE 
WRITE(6, 9560) 

END IF 
C 

IF(JTC.EQ.2) THEN 
C 

C RECALL MATL.NO. AND FIBER DIRECTION ANGLES (FIRST BRANCH) 

C 

MN-MNl (I, J,K) 

A1-ANG1A{I, J,K) 

A2-ANG1B(I, J,K) 

C 

C GET STRESS/ STRAIN MATRIX IN GLOBAL COORDINATES 
C 

CALL GETSS(MM,MN,PROP,SS) 

CALL TRANS2 (A1,A2,T) 

DO 7040 11-1,6 
DO 7040 JJ-1,6 
SUM-0 . 0 

DO 7035 KK-1,6 

7035 SUM-SUM+SS(II,KK) *T(KK, JJ) 

7040 PR(II, JJ)-SUM 
DO 7045 11-1,6 
DO 7045 JJ-1,6 
SUM-0.0 

DO 7044 KK-1,6 

7044 SUM-SUM+T(KK, II)*PR(KK,JJ) 

7045 SSKII, JJ)-SUM 
C 

C RECALL MATL.NO. AND FIBER DIRECTION ANGLES (SECOND BRANCH) 

C 

MN-MN2(I, J,K) 

A1-ANG2A(I, J,K) 

A2-ANG2B(I, J,K) 

C 

C GET STRESS/STRAIN MATRIX IN GLOBAL COORDINATES 
C 

CALL GETSS(MM,MN,PROP,SS) 

CALL TRANS2 (A1,A2,T) 

DO 7050 11-1,6 
DO 7050 JJ-1,6 
SUM-0 . 0 

DO 7049 KK-1,6 

7049 SUM-SUM+SS{II,KK) *T(KK, JJ) 

7050 PR(II, JJ)-SUM 
DO 7055 11-1,6 
DO 7055 JJ-1,6 
SUM-0.0 

DO 7054 KK-1,6 

7054 SUM-SUM+T(KK, II) *PR(KK, JJ) 

7055 SS2(II, JJ)-SUM 
C 

C RECALL FIRST BRANCH VOL. TRACT. AND INTERFACIAL NORMAL ANGLES 
C 

Vl-FVKI, J,K) 

V2-1.0-V1 
A1-AGN1A(I, J,K) 
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A2-AGN1B(I, J,K) 

c 

C GET S.S. MATRICES IN INTEREACIAL, COORDINATES 
C 

CALL TRANS 1 (A1,A2, T; 

DO 7060 11-1,6 
DO 7060 JJ-1,6 
SUM-0.0 
TUM-0 . 0 

DO 7059 KK-1,6 
SUM-SUM-t-SSl (II, KK) ‘KFCK, JJ» 

705 7 TUM-TUM+SS2 (II, KK) *T(KK, JJ) 

RP(II, JJ)-SUM 
7060 PR(II, JJ)-TUM 
DO 7065 11-1,6 
DO 7065 JJ-1,6 
SUM-0.0 
TUM-0 .0 

DO 7064 KK-1,6 
SUM-SUM+T (KK, II) *RP (KK, JJ) 

7064 TUM-TUM+T(KK, II) *PR(KK, JJ) 

SSSKII, JJ)-SUM 

7065 SSS2 (II, JJ)-TUM 
C 

C DO REPLACEMENT MATERIAL ANALYSIS AT FIRST JUNCTION 

C 

CALL GETDD(SSS1,SSS2, V! , ‘.'2, DDl, DD2) 

DO 7075 11-1,6 
DO 7075 JJ-1,6 
DUM-0 .0 

DO 7074 KK-1,6 

7074 DUM-DUM+V1*SSS1 (1I,KK)*DD1{KK,JJ) +V2*SSS2 (II, KK) *DD2 (KK, JJ) 

7075 DD(II, JJ)-DUM 
C 

C GET REPLACEMENT S.S. MATRIX IN GLOBAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 

DO 7079 11-1,6 
DO 7079 JJ-1,6 
SUM-0.0 

DO 7078 KK-1,6 

7078 SUM-SUM+DD(II,KK)*T(KK, JJ) 

7079 PR(II, JJ)-SUM 
DO 7085 11-1,6 
DO 7085 JJ-1,6 
SUM-0.0 

DO 7084 KK-1,6 

7084 SUM-SUM+T(KK, II) *PR(KK, JJ) 

7085 SSKII, JJ)-SUM 
C 

C RECALL MATL.NO. AND FIBER DIRECTION ANGLES (THIRD BRANCH) 

C 

MN-MN3(I, J,K> 

A1-ANG3A(I, J,K) 

A2-ANG3B(I, J,K) 

C 

C GET STRESS/STRAIN MATRIX IN GLOBAL COORDINATES 
C 

CALL GETSS(MM,MN,PROP,SS) 

CALL TRANS2(A1,A2,T) 

DO 7150 11-1,6 
DO 7150 JJ-1,6 
SUM-0.0 

DO 7149 KK-1,6 

7149 SUM-SUM+SS(II,KK) *T(KK, JJ) 

7150 PR(II, JJ)-SUM 
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DO 7155 II-j., 6 
DO 7155 JJ-1,6 
s;w-o . 0 

DO 7.54 KK-1,6 

7154 SOh-SUM+TtKK, II! *PR(KK, .-Ji 

7155 S32 ii^. JJ)-SUM 

C BECADD SECOND BRANCH VO.-.FRACT. ANL INTERFACIAL NORMAL ANGLES 
C 

V2-FV2 (I, J,K) 

A1-AGN2A(I, J,K) 
a2-AGN2B(I, J, K) 

C 

C GET S S. MATRICES IN INTERFACIAL COORDINATES 
C 

CALL TRANSl (Al, A2,TJ 
DO 7160 11-1,6 
DO 7160 JJ-1,6 
;5UM-0 .0 
TUM-0 . 0 

DO 7159 KK-1,6 
SUM-SUM+SSl (II,KK) • T iKK, JJ) 

.59 T5!M-T0M+SS2 (II, KK) *T(KK, JJ) 

HRili, JJ)-SUM 
7160 PP ill, JJ)-TUM 
DC 7165 11-1,6 
DO / 1 6 5 J J— 1,6 
SUM-0.0 
TUM-0 . 0 

DO 7164 KK-1,6 
SUM-SUM+T(KK, II) *RP (KK, JJ) 

7164 TUM-TUM+T (KK, II ) *PR (KK, JJ) 

SSSl (II, JJ)-SUM 

7165 SSS2(II, JJ)-TOM 
C 

C DO REPLACEMENT MATERIAL ANALYSIS AT SECOND JUNCTION 

C 

CALL GETDD(SSS1,SSS2,V;,V2,DD1,0D2) 

C 

C GET AVERAGE STRAINS IN INTERFACIAL COORDINATES 
C 

Ai-AGN2A(I, J,K) 

A2-AGN2B(I, J,K) 

CALL TRANS2 (A1,A2,T) 

DO 7250 11-1,6 
DOM-0 .0 

DO 7240 JJ-1,6 
7240 DUM-DUM+T(II, JJ)*TS(JJ) 

7250 TSS(II)-DUM 
C 

C GET CONSTITUENT STRAINS IN INTERFACIAL COORDIDATES 
C 

DO 7311 11-1,6 
TSl(II)-0.0 
TS2(II)-0.0 
DO 7311 JJ-1,6 

TSl (II)-TS1(II)+DD1(II, JJ)*TSS(JJ! 

7311 TS2(I1)-TS2(II)-^DD2(I1, JJ)*TSS{JJ) 

C 

C GET CONSTITUENT STRAINS IN GLOBAL COORDINATES 
C 

CALL TRANSl (A1,A2,T) 

DO 7320 11-1,6 
DUM-0 .0 
TUM-0 . 0 
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DO 7315 JJ-1, 6 
DUM-DUM+T (II, JJ) •T':' - . - 
73.5 TUM-TUM+T ( 1 1 , JJ> ■ 

TS (II)-DUM 
/320 TSS (II) -TUM 

A1-ANG3A (I, J, K) 

A2-ANG3B(I, J, K) 

CALL TRANS2 (Al, A2. T; 

DO 7330 11-1,6 
DUM-0 . 0 

DO 7325 JJ-1, 6 
7i_5 DUM-DUM+T(II,JJ)*TSS 
jO TS2(II)-DUM 

C GET COWSTITOENT STRESSES IN MAVGPlHi, 30JiiGGNATSS 

c 

MN-MN3 ( 1, J, K) 

CALL GETSS (MM,MN, fROP, '.S, 

CALL GETMS (MM,MN, ?ROP,'TS2,£»^f'E, 

DO 7350 II-l, 6 
ST2 (II)-O.O 
DO 7350 JJ-1, 6 

-,S0 ST2 (II) -ST2 (II) +SS (II, . J) '-'S2 ( JJ 
WRITE i6, 9560) 

WRITE (6, 9690) I,J,K 
WRITE (6, 9692) MN 

WRITE(6, 9694) st: (1) , ST2 i : i ,S-2 ( ,i; 

WRITE (6, 9696) ST2{4),f: ■ 5; . .STJ {■,( 

WRITE (6, 9698) SAFE 
WRITE l6, 9560) 

C 

t recall first branch InTEV. -."Io ncrma^ angles 

A1-AGN1A(I, J,K) 

A2-AGN1B(I, J,K) 

: GET AVERAGE STRAINS IN INTEHEAC 1.AL COORDINATES 

CALL TRANS2 (Al, A2,T) 

DO 7370 II-l, 6 
DUM-0 . 0 

DO 7360 JJ-1, 6 
7360 DUM-DUM+T (II, JJ) *TS (JJ) 

7370 TSS(II)-DUM 

r RECALL FIRST BRANCH MATL. NUMBER A?.'D FIBER ANGLES 

c; 

MN-MNl (I, J,K) 

A1-ANG1A(I, J, K) 

A2-ANG1B(I, J, K) 

c 

C GET FIRST BRANCH S.S. MATRIX IN GLOBAL COORDINATES 

C 

CALL GETSS(MM,MN. ?ROP,S": 

CALL TRANS2 (Al, A2,T, 

DO 7400 11-1,6 
DO 7400 JJ-l, 6 
SUM-0.0 

DO 7390 KK-1,6 

7390 SUM-SUM+SS(II,KK) *T(KK, JJ) 

7400 PR(II, JJ)-SUM 
DO 7450 11-1,6 
DO 7450 JJ-1, 6 
SUM-0.0 

DO 7440 KK-1,6 

7440 SUM-SUM+T(KK, II) *PR(KK, JJ) 
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7450 SSI (II, JJ) -SUM 

C RECMiL SECOND BRANCH MAIL- NUMBER AND FIBER ANGLES 
C 

HN-MN2(I, J,K) 

A1-ANG2A(I, J,K) 

A2-ANG2B(I, J,K) 

C GET SECOND BRANCH S.S. MATRIX IN GLOBAL C(X)RDINATES 
C 

CALL GETSS(MM,MN,PROP,SS) 

CALL TRANS2 (A1,A2,T) 

DO 7500 11-1,6 
DO 7500 JJ-1,6 
SUM-0.0 

DO 7490 KK-1,6 

7490 SUM-SUM+SS (II,KK) *T(KK, JJ) 

7500 PR(II, JJ)-SUM 
DO 7550 11-1,6 
DO 7550 JJ-1,6 
SUM-0 . 0 

DO 7540 KK-1,6 

7540 SUM-SUM+T(KK, II)*PR(KK, JJ) 

7550 SS2(II,JJ)-SUM 

C RECALL FIRST BRANCH VOL. FRACT. AND INTER. NORMAL ANGLES 
C 

Vl-FVld, J,K> 

V2-1.0-V1 
A1-AGN1A(I, J,K) 

A2-AGN1B(I, J,K) 

C GET MATL S.S. MATRICES IN INTERFACIAL COORDINATES 
C 

CALL TRANS1(A1,A2,T) 

DO 7600 11-1,6 
DO 7600 JJ-1,6 
SUM-0 . 0 
TUM-0.0 

DO 7590 KK-1,6 

SUM-SUM+SS1(II,KK)*T(KK,JJ) 

7590 TUM-TUM+SS2 (II,KK) *T(KK, JJ) 

RP(II,JJ)-SUM 
7600 PR(II, JJ)-TUM 
DO 7650 11-1,6 
DO 7650 JJ-1,6 
SUM-0 . 0 
TUM-0 . 0 

DO 7640 KK-1,6 
SUM-SUM+T (KK, II ) *RP (KK, JJ) 

7640 TUM-TUM+T(KK, II)*PR(KK, JJ) 

SSS1(II,JJ) -SUM 
7650 SSS2(II, JJ)-TUM 
C 

C DO REPLACEMENT MATERIAL ANALYSIS 
C 

CALL GETDD(SSS1,SSS2,V1,V2,DD1,DD2) 

Q 

C GET CONSTITUENT STRAINS IN INTERFACIAL COORDINATES 
C 

DO 7711 11-1,6 
TSl(II)-0.0 
TS2(II)-0.0 
DO 7711 JJ-1,6 

TS1(II)-TS1 (II)+DD1 (II, JJ) *TSS(JJ) 

7711 TS2 (II) -TS2 (II) +DD2 (II, JJ) *TSS ( JJ) 
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c 

C GET CONSTITUENT STRAINS IN GLOBAL COORDINATES 
C 

DO 7720 11-1,6 

DOM-0.0 

TUM-0.0 

DO 7715 JJ-1,6 
DOM-DUM+T <11, JJ) *TS1 { JJ) 

7715 TUM-TUM+T(II, JJ) *TS2(JJ) 

TS(II)-DUM 
7720 TSS(I1)-TUM 
C 

C RECALL FIRST BRANCH FIBER ANGLES 
C 

A1-ANG1A(I, J,K) 

A2-ANG1B(I, J,K) 

C 

C GET FIRST BRANCH STRAINS IN MATERIAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 

DO 7730 11-1,6 
DOM-0.0 

DO 7725 JJ-1,6 
7725 DUM-DUM+T(II, JJ)*TS(JJ) 

7730 TS1(II)-DUM 
C 

C RECALL SECOND BRANCH FIBER ANGLES 
C 

A1-ANG2A(I, J,K) 

A2-ANG2B(I, J,K) 

C 

C GET SECOND BRANCH STRAINS IN MATERIAL COORDINATES 
C 

CALL TRANS2 (A1,A2,T) 

DO 7740 11-1,6 
DUM-0.0 

DO 7735 JJ-1,6 

7735 DOM-DUM+T (II, JJ)*TSS(JJ) 

7740 TS2(II)-DOM 
C 

C GET FIRST BRANCH STRESSES 
C 

MN-MNl (I, J,K) 

CALL GETSS(MM,MN,PROP,SS) 

CALL GETMS(MM,MN, PROP, TSl, SAFE) 

DO 7750 11-1,6 
STl(II)-0.0 
DO 7750 JJ-1,6 

7750 STl (II) -STl (II) +SS (II, JJ) *TS1 (JJ) 

WRITE (6, 9560) 

WRITE(6,9690) I,J,K 
WRITE (6, 9692) MN 

WRITE (6, 9694) STl (1) , STl (2) , STl (3) 

WRITE (6, 9696) STl (4) , STl (5) , STl (6) 

WRITE (6, 9698) SAFE 
WRITE (6, 9560) 

C 

C GET SECOND BRANCH STRESSES 
C 

MN-MN2(I, J,K) 

CALL GETSS(MM,MN,PROP,SS) 

CALL GETMS (MM, MN, PROP, TS2, SAFE) 

DO 7760 11-1,6 
ST2(Il)-0.0 
DO 7760 JJ-1,6 

7760 ST2 (II) -ST2 (II) +SS(II, JJ) *TS2 (JJ) 
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WRITE (6, 9560) 

WRITE (6, 9690) I,J,K 
WRITE<6,9692) MW 

WRITE (6, 9694) ST2 (1) ,ST2 (2) ,ST2 (3) 
WRITE (6, 9696) ST2 (4 ) , ST2 (5) ,ST2 (6) 
WRITE (6, 9698) SAFE 
WRITE (6, 9560) 

END IF 

8000 CONTINUE 
C 

9000 FORMAT (F12. 5) 

9010 FORMAT (F12. 2) 

9015 FORMAT (E16. 6) 

9020 FORMAT (6F12. 2) 

9025 FORMAT (10X,4E14. 4) 

9030 FORMAT (15) 


C9031 FORMAT (315) 
9060 FORMAT (IH , 
9080 FORMAT (IH , 
9090 FORMAT (IH , 
9095 FORMAT (IH , 
9100 FORMAT (IH , 
9120 FORMAT (IH , 
9130 FORMAT (IH , 
9140 FORMAT (IH , 
9150 FORMAT (IH , 
9160 FORMAT (IH , 
9170 FORMAT (IH , 
9175 FORMAT (IH , 

9177 FORMAT (IH , 

9178 FORMAT (IH , 

9179 FORMAT (IH , 

9180 FORMAT (IH , 
9190 FORMAT (IH , 
9320 FORMAT (IH , 
9440 FORMAT (IH , 
9450 FORMAT (IH , 
9455 FORMAT (IH , 
9460 FORMAT (IH , 
9480 FORMAT (IH , 
9485 FORMAT (IH , 
9490 FORMAT (IH , 
9495 FORMAT (IH , 
9500 FORMAT (IH , 
9510 FORMAT (IH , 
9520 FORMAT (IH « 
9525 FORMAT (IH , 
9530 FORMAT (IH , 
9540 FORMAT (IH , 
9550 FORMAT (IH , 


'INPUT NUMBER OF JUNCTIONS AT LOCATION' , 314) 
'INPUT NO. SUBCELLS (X DIR.) IN UNIT CELL') 

'INPUT NO. SUBCELLS (Y DIR.) IN UNIT CELL') 

'INPUT NO. SUBCELLS (Z DIR.) IN UNIT CELL') 

'INPUT NO. COMPOSITE MATERIALS NEEDED, NM') 

' INPUT E IN FIBER DIRECTION' ) 

' INPUT E NORMAL TO FIBER DIRECTION' ) 

' INPUT MAJOR POISSONS RATIO IN LT PLANE' ) 
'INPUT POISSONS RATIO IN TT PLANE') 

'INPUT SHEAR MODULUS G IN LT PLANE') 

'INPUT SHEAR MODULUS G IN TT PLANE') 

'INPUT LONG. TENSION ALLOWABLE') 

'INPUT LONG COMPRESSION ALLOWABLE') 

'INPUT TRANS. TENSION ALLOWABLE') 

'INPUT TRANS. COMPRESSION ALLOWABLE') 

'SELECT A MATERIAL NUMBER FROM ONE TO TEN') 
'MATERIAL PROPERTY DATA ECHO') 

'SPECIFY THE CURRENT MATL. ID. NO.') 

'INPUT SIDE LENGTH OF UNIT CELL IN X DIR.') 

'INPUT SIDE LENGTH OF UNIT CELL IN Y DIR.') 

'INPUT SIDE LENGTH OF UNIT CELL IN Z DIR.') 

'INPUT DIST.(%) ORIGIN TO UNIT CELL NODE', 13) 
'INPUT 1ST FIBER SPHERICAL ANGLE') 

'INPUT 1ST INTERFACIAL NORMAL ANGLE') 

'INPUT 2ND FIBER SPHERICAL ANGLE') 

'INPUT 2ND INTERFACIAL NORMAL ANGLE') 
'EX,EY,EZ - ',3F12.2) 

'GYZ,GXZ,GXY - ',3F12.2) 

'MUYZ,MUXZ,MUXY - ',3F12.4) 

'MUZY,MUZX,MUYX - ',3F12.4) 

'NUYZ,X ; NUYZ,Y ; NUYZ, Z - ',3F12.4) 

'NUXZ,X ; NUXZ,Y ; NUXZ, Z - ',3F12.4) 

'NUXY,X ; NUXY,Y ; NUXY, Z - ',3F12.4) 


9560 FORMAT (IH ) 

9600 FORMAT (IH , 13X, ' ELASTIC CONSTANTS OF THE COMPOSITE ') 

9610 FORMATdH , 13X, ' INPUT APPLIED STRESSES IN X,Y,Z COORDINATES') 
9620 FORMATdH ,5X,' INPUT X NORMAL STRESS ') 

9630 FORMATdH , 5X, ' INPUT Y NORMAL STRESS ') 

9640 FORMATdH ,5X,' INPUT Z NORMAL STRESS ') 

9650 FORMATdH ,5X,' INPUT YZ SHEAR STRESS ') 

9660 FORMATdH , 5X, ' INPUT XZ SHEAR STRESS ') 

9670 FORMATdH ,5X,' INPUT XY SHEAR STRESS ') 

9690 FORMATdH , 5X, ' STRESSES IN ELEMENT NO. ',313) 

9692 FORMATdH , 'MATERIAL NO. ',13) 

9694 FORMATdH ,' NORMAL 1,2,3 ',3F14.2) 

9696 FORMATdH , 'SHEAR 23,13,12 ',3F14.2) 

9698 FORMATdH , 'MINIMUM MARGIN OF SAFETY IS ',F12.4) 

9700 FORMATdH ,' INPUT MATL. NO. 1 AT ',314) 

9710 FORMATdH ,' INPUT MATL. NO. 2 AT ',314) 
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ono 


9715 FORMAT (IH 
9720 FORMAT (IH 
9895 FORMAT (IH 

9898 FORMAT (IH 

9899 FORMAT (IH 


INPUT MATL. NO. 

, • INPUT 1ST MATL. 

((',F7.1,'))') 
,» ((',F4.1,'))') 

' ((M3,' ))') 


3 AT ',314) 

VOLUME FRACTION') 


9999 END 

MATRIX INVERSION BY PARTITIONING OF A 6X6 MATRIX 


SUBROUTINE INV(A) 

DIMENSION A(6,6),B(4,4),X(4,4),Y(4,4),Z(4,4) 
D-A(l, 1) * (A(2,2)*A(3, 3) -A(2, 3) *A(3, 2) ) 

1 -A(2,l) * (A(1,2)*A(3,3)-A(1,3)*A{3,2) ) 

2 +A(3,1)MA(1,2)*A(2,3)-A(2,2)*A(1,3)) 
IF(D.EQ.O.O) GOTO 700 

B(1,D- (A(2,2)*A(3,3)-A(2,3)*A 3,2 /D 

B(l, 2) — (A(l, 2) *A(3, 3)-A(l, 3) *A(3,2) ) /D 
B(1,3)- (A(1,2)*A(2,3)-A(1,3)*A(2,2))/D 
B(2,D — (A(2,1)*A(3,3)-A(2,3)*A(3,1))/D 
B(2,2)- (A(1,1)*A(3,3)-A(1,3)*A(3,1))/D 
B(2,3) — (A(1,1)*A(2,3)-A(1,3)*A(2,1))/D 
B(3,1>- (A(2,l)*A(3,2)-A(2,2)*A(3,l))/D 
B(3,2) — (A(l, 1) *A(3,2)-A(1,2) *A(3,1))/D 

B(3,3)- (A(1,1)*A(2,2)-A(1,2)*A(2,1))/D 

DO 100 1-1,3 

DO 100 J-1,3 

X(I, J)-0.0 

Y(I, J)-0.0 

100 Z(I,J)-A(I+3, J+3) 

DO 200 1-1,3 
DO 200 J-1, 3 
DO 200 K-1,3 

X(I,J)-X(I,J)+B(I,K)*A(K,J+3) 

200 Y(I, J)-Y(I, J)+A(I+3,K)*B(K, J) 

DO 300 1-1,3 
DO 300 J-1,3 
DO 300 K-1,3 

300 Z(I,J)-Z(I,J)-T(I,K)*A(K,J+3) 

D-Z(1,1)*(Z(2,2)*Z(3,3)-2(2,3)*Z(3,2)) 

1 -Z(2,1)*(Z(1,2)*Z(3,3)-Z(1,3)*Z(3,2)) 

2 +Z(3,D* (Z(1,2)*Z(2,3)-Z(2,2)*Z(1,3) ) 
IF(D.EQ.O.O) GOTO 700 


DO 400 1-1,6 
DO 400 J-1, 6 
400 A(I,J)-0.0 

A (4, 4)- (Z(2,2)*Z(3,3)-Z(2,3)*Z(3,2))/D 
A(4,5) — (Z(1,2)*Z(3,3)-Z(1,3)*Z(3,2))/D 
A (4, 6)- (Z(1,2)*Z(2,3)-Z(1,3)*Z(2,2))/D 
A(5,4) — (Z(2,1)*Z(3,3)-Z(2,3)*Z(3,1))/D 
A (5, 5)- (Z(1,1)*Z(3,3)-Z(1,3)*Z(3,1))/D 
A(5,6) — (Z(1,1)*Z(2,3)-Z(1,3)*Z(2,1))/D 
A(6,4)- (Z(2,1)*Z(3,2)-Z(2,2)*Z(3,1))/D 
A (6, 5) — (Z (1,1) *Z(3,2)-Z(1,2)*Z(3,1))/D 
A(6,6)- (Z(1,1)*Z(2,2)-Z(1,2)*Z(2,1))/D 
DO 500 1-1,3 
DO 500 J-1,3 
A(I,J)-B(I, J) 

DO 500 K-1,3 

A(I,J+3)-A(I, J+3)-X(I,K) *A(K+3,J+3) 

500 A (1+3 , J) -A ( 1+3, J) -A ( 1+3 , K+3) *Y (K, J) 

DO 600 1-1,3 
DO 600 J-1,3 
DO 600 K-1,3 

600 A(I,J)-A(I, J)-A(I,K+3)*Y(K,J) 

GO TO 2000 ^ . 

700 D- A(1,1)*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A{4,3)) 
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1 -A(2,3)*(A(3,2)*A(4,4)-A(3,4)*A(4,2)) 

2 +A(2,4)MA(3,2)*A(4,3)-A(3,3)*A(4,2))) 
D-D-A(1,2)*{A(2,1)*(A{3,3)*A(4,4)-A{3,4)*A<4,3)) 

1 -A(2,3)*(A(3,1)*A{4,4)-A(3,4)*A(4,1)) 

2 +A(2,4)MAP,1)*A{4,3)-A(3,3)*A(4,1))) 
D-D+A(1,3)*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2)) 

1 -A(2,2)*(A(3,1)*A(4,4)-A(3,4)*A(4,1>) 

2 ♦A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1))) 
D-D-A(1,4)MA(2,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2)> 

1 -A(2,2)*(A(3,1)*A(4,3)-A(3,3)*A(4,D) 

2 +A(2,3)*(A{3,1)*A(4,2)-A(3,2)*A(4,1))> 
C WRITE (6, 9010) D 

C9010 FORMAT (F12. 2) 

IF(D.EQ.O.O) D-1.0 




1 

2 

I 

1 

2 

I 

1 

2 


B(l,2)- 


B(l,3)- 


B(l,4)“ 


1 

2 

I 

1 

2 

I 

1 

2 

I 

1 

2 

1 

1 

2 

1 

1 

2 

I 

1 

2 

1 

1 

2 


B(2,D- 


B(2,2)- 


B(2,3)- 


B(2,4)' 


B(3,D- 


B(3,2)- 


B(3,3)- 


B(3,4)- 


1 

2 

I 

1 

2 

! 

1 

2 

I 

1 

2 

I 

1 

2 


B(4,D- 


B(4,2)- 


B(4,3)- 


B(4,4)- 


1090 


DO 

DO 


1090 

1090 


B(l, J)-B(I, J)/D 
DO 1100 1-1,4 


-A(2,3)MA(3,2)*A 
+A(2,4)*(A(3,2)*A 
-A(1,2)*(A(3,3)*A 
+A(1,3)*(A{3,2)*A 
-A(1,4)*(A(3,2)*A 
+A(1,2)*(A{2,3)*A 
-A(1,3)*(A(2,2)*A 
+A(1,4)*(A{2,2)*A 
-A(1,2)*(A(2,3)*A 
+A{1,3)MA(2,2)*A 
-A{1,4)*(A(2,2)*A 
-A(2,1)*(A(3,3)*A 
+A(2,3)*(A(3,1)*A 
-A(2,4)*(A(3,1)*A 
+A(1,1)*(A(3,3)*A 
-A{1,3)*(A(3,1)*A 
+A(1,4)*(A(3,1)*A 
-A(1,1)*(A<2,3)*A 
+A(1,3)*(A(2,1)*A 
-A(1,4)*(A<2,1)*A 
+A(1,1)*(A(2,3)*A 
-A(1,3)*(A(2,1)*A 
+A(1,4)MA(2,1)*A 
+A(2,1>*(A(3,2)*A 
-A(2,2)*(A(3,1)*A 
+A(2,4)*(A(3,1)*A 
-A(1,1)MA(3,2)*A 
+A(1,2)MA(3,1)*A 
-A(1,4)*(A(3,1)*A 
+A(1,1)*<A(2,2)*A 
-A(1,2)*(A(2,1)*A 
+A(1,4)*(A(2,1)*A 
-A(1,1)*(A(2,2)*A 
+A(1,2)*(A(2,1)*A 
-A(1,4)*(A(2,1)*A 
-A<2,1)*(A(3,2)*A 
+A(2,2)MA(3,1)*A 
-A{2,3)MA(3,1)*A 
+A(1,1)MA(3,2)*A 
-A(1,2)*(A(3,1)*A 
+A(1,3)MA(3,1)*A 
-A(1,1)*(A(2,2)*A 
+A(1,2)MA(2,1)*A 
-A(1,3)*<A(2,1)*A 
+A(1,1)*(A(2,2)*A 
-A(1,2)MA(2,1)*A 
+A(1,3)*(A(2,1)*A 
1-1,4 
J-1,4 


(4,4)-A(3,4 

*A{4,3)> 

(4,4)-A(3,4 

*A(4,2)) 

(4,3)-A(3,3 

*A(4,2)) 

(4,4)-A(3,4 

*A{4,3)> 

(4,4)-A(3,4 

*A(4,2)) 

(4,3)-A(3,3 

*A(4,2)) 

(4,4)-A(2,4 

*A(4,3)) 

(4,4)-A<2,4 

*A(4,2)) 

(4,3)-A<2,3 

*A(4,2)) 

(3,4>-A(2,4 

*A(3,3)) 

(3,4)-A(2,4 

*A(3,2)) 

(3,3)-A(2,3 

*A{3,2)) 

(4,4)-A(3,4 

*A(4,3)) 

(4,4)-A(3,4 

*A(4,D) 

(4,3)-A(3,3 

*A{4,D) 

(4,4)-A(3,4 

*A(4,3)) 

(4,4)-A{3,4 

*A<4,D) 

(4,3)-A(3,3 

*A<4,D) 

(4,4)-A(2,4 

*A{4,3)) 

(4,4)-A(2,4 

*A(4,D) 

(4,3)-A(2,3 

*A(4,D) 

(3,4)-A(2,4 

*A(3,3)) 

(3,4)-A(2,4 

*A(3,D) 

(3,3)-A(2,3 

*A{3,D) 

(4,4)-A<3,4 

*A(4,2)> 

(4,4)-A(3,4 

*A(4,D) 

(4,2)-A{3,2 

*A(4,D) 

(4,4)-A{3,4 

*A(4,2)) 

(4,4)-A(3,4 

*A(4,D) 

(4,2)-A(3,2 

*A(4,D) 

(4,4)-A{2,4 

*A{4,2)) 

(4,4)-A(2,4 

*A(4,D) 

(4,2)-A(2,2 

*A(4,D) 

(3,4)-A(2,4 

*A{3,2)) 

(3,4)-A(2,4 

*A{3,1)> 

(3,2)-A(2,2 

*A(3,D) 

,(4,3)-A(3,3 

*A(4,2)) 

,(4,3)-A<3,3 

*A(4,D) 

,{4,2)-A<3,2 

*A(4,D) 

,{4,3)-A(3,3 

*A{4,2)) 

,(4,3)-A{3,3 

*A(4,D) 

,(4,2)-A(3,2 

*A(4,D) 

,(4,3)-A(2,3 

*A{4,2)) 

.(4,3)-A(2,3 

*A{4,1>) 

>(4,2)-A(2,2 

*A(4,D) 

.{3,3)-A(2,3 

*A(3,2)) 

i(3,3)-A(2,3 

*A(3,D) 

l(3,2)-A(2,2 

*A(3,D) 
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DO HOO J-1,4 
X(I, J)-0.0 

1100 Y{I,J)-0.0 

DO 1110 1-1,2 
DO 1110 J-1,2 

1110 Z(I,J)-A(I+4, J+4) 

DO 1200 1-1,4 
DO 1200 J-1,2 
DO 1200 K-1,4 

1200 X(I,J)-X(I, J)+B(I,K)*A(K, J+4) 

DO 1210 1-1,2 
DO 1210 J-1,4 
DO 1210 K-1,4 

1210 Y(I, J)-Y<1, J)+A(I+4,K)*B(K, J) 

DO 1300 1-1,2 
DO 1300 J-1,2 
DO 1300 K-1,4 

1300 Z(I, J)-Z(I, J)-Y(I,K)*A(K, J+4) 
D-Z(1,1)*Z(2,2)-Z(1,2)*Z(2,1) 

XF(D.EQ.O.O) D-1.0 
DO 1400 1-1,6 
DO 1400 J-l,€ 

1400 A(I,J)-0.0 

A(5,5)-Z(2,2)/D 
A<5,6)— Z(l,2)/D 
A(6,5)— Z(2,l)/D 
A(6,6)-Z(l,l)/D 
DO 1450 1-1,4 
DO 1450 J-1,4 

1450 A(I, J)-B<I, J) 

DO 1500 1-1,4 
DO 1500 J-1,2 
DO 1500 K-1,2 

1500 A(I, J+4)-A(I, J+4)-X(I,K)*A(K+4, J+4) 

DO 1550 1-1,2 
DO 1550 J-1,4 
DO 1550 K-1,2 

1550 A(I+4, J)-A(I+4, J)-A(I+4,K+4)*Y(K, J) 

DO 1600 1-1,4 
DO 1600 J-1,4 
DO 1600 K-1,2 

1600 A(I, J)-A(I, J)-A{I,K+4)*Y(K, J) 

2000 CONTINUE 
RETURN 
END 

C 

SUBROUTINE MATINV (A, NMAX,N,B, MAX, M, DETERM) 

C IMPLICIT REAL*8 (A-H,0-Z) 

C**** STANDARD MATRIX INVERSION SUBPROGRAM 

C**** 

DIMENSION A(NMAX,NMAX),B<NMAX,MAX) 

DIMENSION IPIVOT (300) , INDEX (300,2), PIVOT (300) 

INITIALIZATION 

10 DETERM - 1.0 
15 DO 20 J-1,N 
20 IPIVOT (J) - 0 
30 DO 550 1-1, N 

SEARCH FOR THE PIVOT ELEMENT 

40 AMAX - 0.0 

45 DO 105 J-1,N 

50 IF (IPIVOT (J) .EQ. 1) GOTO 105 
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60 

DO 100 K-1,N 

70 

IF (IPIVOT 

80 

IF (ABS (J 

85 

IROW - 

90 

ICOLUM 

95 

AMAX - 

100 

CONTINUE 

105 

CONTINUE 

110 

IPIVOT (ICOLUM) - 


1)00,100,740 


C INTERCHANGE ROWS TO PUT ELEMENT ON DIAGONAL 

^ 130 IF (IRON .EQ. ICOLUM) GOTO 260 
140 DETERM - -DETERM 

150 DO 200 L-1,N 

160 SWAP - A(IROW,L) 

170 A(IROW,L) - A (ICOLUM, L) 

200 A (ICOLUM, L) - SWAP 

205 IF (M .LE. 0) GOTO 260 

210 DO 250 L-1,M 

220 SWAP - B(IROW,L) 

230 B(lROW,L) - B( ICOLUM, L) 

250 B (ICOLUM, L) - SWAP 

260 INDEX (1,1) - IROW 

270 INDEX (I, 2) - ICOLUM 

310 PIVOT (I) - A (ICOLUM, ICOLUM) 

320 DETERM - DETERM*PIVOT (I) 


DIVIDE PIVOT BY PIVOT ELEMENT 

330 A (ICOLUM, ICOLUM) - 1,0 
340 DO 350 L-1,N 

350 A ( ICOLUM, L) - A (ICOLUM, L) /PIVOT (I) 
355 IF (M ,LE. 0) GOTO 380 
360 DO 370 L-1,M 

370 B( ICOLUM, L> - B (ICOLUM, L) /PIVOT (I) 


C 

C REDUCE NON-PIVOT ROWS 
C 

380 DO 550 L1-1,N 

390 IF (LI ,EQ. ICOLUM) GOTO 550 

400 T - A (LI, ICOLUM) 

420 A (LI, ICOLUM) - 0.0 

430 DO 450 L-1,N 

450 A(L1,L) - A(L1,L) - A(ICOLUM,L) *T 

455 IF (M ,LE. 0) GOTO 550 

460 DO 500 L-1,M 

500 B(L1,L) - B(L1,L) - B (ICOLUM, L) *T 

550 CONTINUE 


C 

C INTERCHANGE COLUMNS 
C 

600 DO 710 I-1,N 
610 L - N + 1 - I 

620 IF (INDEX(L,1) .EQ. INDEX(L,2)) GOTO 710 
630 JROW - INDEX (L,l) 

640 JCOLUM - INDEX(L,2) 

650 DO 705 K-1,N 

660 SWAP - A (K, JROW) 

670 A (K, JROW) - A (K, JCOLUM) 

700 A (K, JCOLUM) - SWAP 

705 CONTINUE 

710 CONTINUE 
740 RETURN 
END 
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SUBROUTINE TRANS 1 (Al» A2,T) 

C GET THE STRAIN TRANSFORMATION MATRIX (T) 
C 

DIMENSION T(6,6) 

51- SIND (Al) 

52- SIND(A2) 

Cl-COSD (Al) 

C2-COSD(A2) 

S1S-S1*S1 

S2S-S2*S2 

C1S-C1*C1 

C2S-C2*C2 

SC1-S1*C1 

SC2-S2*C2 

T(1,1)-C1S*C2S 

T(2,1)-S1S*C2S 

T(3,1)-S2S 

T(4,1)-S1*SC2*2.0 

T(5,1)-C1*SC2*2.0 

T(6,1)-SC1*C2S*2.0 

T(1,2)-S1S 

T(2,2)-C1S 

T(3,2)-0.0 

T(4,2)-0.0 

T(5,2)-0.0 

T(6,2)— SC1*2.0 

T(1,3)-C1S*S2S 

T(2,3)-S1S*S2S 

T(3,3)-C2S 

T(4,3)— S1*SC2*2.0 

T(5,3)— C1*SC2*2.0 

T(6,3)-SC1*S2S*2.0 

T(1,4)-SC1*S2 

T(2,4)— T{1,4) 

T(3,4)-0.0 
T{4,4)-C1*C2 
T(5,4)— S1*C2 
T(6,4) — (C1S-S1S)*S2 
Ta,5)— C1S*SC2 
T(2,5)— S1S*SC2 
T(3,5)-SC2 
T(4,5)-S1*(C2S-S2S) 
T(5,5)-C1*(C2S-S2S) 

T(6,5)— 2.0*SC1*SC2 
T(l,6)— SC1*C2 
T(2,6)-SC1*C2 
T(3,6)-0.0 
T(4,6)-C1*S2 
T(5,6)— S1*S2 
T(6,6)-(C1S-S1S)*C2 
RETURN 

END 


SUBROUTINE TRANS2 (Al, A2,T) 

GET THE INVERSE STRAIN TRANSFORMATION MATRIX (T) 

DIMENSION T(6,6) 

51- SIND(Al) 

52- SIND(A2) 

Cl-COSD (Al) 

C2-COSD(A2) 

S1S-S1*S1 

S2S-S2*S2 

C1S-C1*C1 
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C2S-C2*C2 

SC1-S1*C1 

SC2-S2*C2 

T(1,1)-C1S*C2S 

T(1,2)-S1S*C2S 

T(1,3)-S2S 

T(1,4)-S1*SC2 

T(1,5)-C1*SC2 

T(1,6)-SC1*C2S 

T(2,1)-S1S 

T(2,2)-C1S 

T(2,3)-0.0 

T(2,4)-0.0 

T(2,5)-0.0 

T{2,6)— SCI 

T(3,1)-C1S*S2S 

T(3,2)-S1S*S2S 

T{3,3)-C2S 

T(3,4)— S1*SC2 

T(3,5)— C1*SC2 

T(3,6)-SC1*S2S 

T(4,1)-2.0*SC1*S2 

T(4,2)— T(4,l) 

T(4,3)-0.0 
T(4,4)-C1*C2 
T(4,5)— S1*C2 
T(4,6) — (C1S-S1S)*S2 
T(5,D— 2.0*C1S*SC2 
T(5,2)— 2.0*S1S*SC2 
T(5,3)-2.0*SC2 
T(5,4)-S1*(C2S-S2S) 
T(5,5)-C1*(C2S-S2S) 

T(5,6)— 2.0*SC1*SC2 

T(6, 1)— 2 .0*SC1*C2 

T(6,2)-2.0*SC1*C2 

T(6,3)-0.0 

T(6,4)-C1*S2 

T(6,5)— S1*S2 

T(6,6)-(C1S-S1S)*C2 

RETURN 

END 

C 

SUBROUTINE GETSS (MM,MN, PROP, SS) 

C 

DIMENSION SS( 6, 6), PROP (MM, 10) 

DO 1113 11-1,6 
DO 1113 JJ-1,6 
1113 SS(II,JJ)-0.0 
El-PROP(MN,l) 

E2-PROP (MN,2) 

Ul-PROP (MN, 3) 

U2-PROP(MN,4) 

Gl-PROP <MN,5) 

G2-PROP(MN,6) 

R-U1*U1*E2/E1 

D-(1.0+U2) * (1.0-U2-2.0*R) 

SS(1,1)-E1*(1.0-U2*U2)/D 

SS(1,2)-E2*U1*(1.0+U2)/D 

SS(1,3)-SS{1,2> 

SS(2,1)-SS(1,2) 

SS(2,2)-E2*(1.0-R)/D 

SS(2,3)-E2*(U2+R)/D 

SS(3,1)-SS(1,3) 

SS(3,2)-SS(2,3) 

SS(3,3)-SS(2,2) 

SS(4,4)-G2 
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SS(5,5)-G1 

SS(6,6)-G1 

RETURN 

END 


C 


SUBROUTINE GETB(AA, BB,CC,X, 


DIMENSION BM(6,24) 

DO 100 I-l»6 
DO 100 J-1,24 

100 125/AA) * {1 . 0-Y/BB) * ( 1 . 

Yl— (0 . 125/BB) * (1 ,0-X/AA) * (1. 0-Z/CC) 
21— (0 . 125/CC) * (1 .0-X/AA) * <1 
X2— (0 . 125/AA) * (1 .0-Y/BB) * (1.0+Z/CC) 
Y2“- (0 . 125/BB) * (1.0-X/AA)*(1.0+Z/CC) 
Z2-+ (0 . 125/CC) * (1 .0-X/AA) * (1.0-Y/BB) 
X3— (0 . 125/AA) * (1 . 0+Y/BB) * (1 . 0-Z/CC) 
Y3-+ (0 . 125/BB) * (1 .0-X/AA) * (1 .0-Z/CC) 

23— (0. 125/CC) *(1.0-X/AA)*(1.0+Y/BB) 

X4— (0 . 125/AA) * (1 .0+Y/BB) * (1 . 0+Z/CC) 
Y4-+ (0. 125/BB) *{1. 0-X/AA) * (1.0+Z/CC) 

24- + (0 . 125/CC) * (1 . 0-X/AA) * (1 .0+Y/BB) 
X5-+(0. 125/AA) * (1.0-Y/BB) * (1. 0-Z/CC) 
Y5— (0. 125/BB) *(1.0+X/AA)* (1. 0-Z/CC) 

25— (0 . 125/CC) * (1 .0+X/AA) * (1 ,0-Y/BB) 
X6-+ (0 . 125/AA) * (1.0-Y/BB) * (1.0+Z/CC) 
Y6—(0. 125/BB) * (1. 0+X/AA) *(1.0+2/CC) 

26- + (0 . 125/CC) * (1 . 0+X/AA) * (1 . 0-Y/BB) 
X7-+ { 0 . 125/AA) * (1 . 0+Y/BB) * ( 1 . 0-Z/CC) 
Y7-+ (0 . 125/BB) * (1 .0+X/AA) * (1 .0-Z/CC) 

27— (0 . 125/CC) * (1 .0+X/AA) * (1. 0+Y/BB) 
X8-+ (0 . 125/AA) * (1 .0+Y/BB) * (1.0+Z/CC) 
Y8-+ (0 . 125/BB) * (1 .0+X/AA) * (1.0+Z/CC) 

28- + (0 . 125/CC) * (1 . 0+X/AA) * (1 .0+Y/BB) 
BM(1,1)-X1 

BM(2,2)-Y1 

BM(3,3)-Z1 

BM(1,4)-X2 

BM(2, 5)-Y2 

BM(3,6)-Z2 

BM(1,7)-X3 

BM(2,8)-Y3 

BM(3,9)-Z3 

BM(1,10)-X4 

BM(2,11)-Y4 

BM(3,12)-Z4 

BM(1,13)-X5 

BM(2,14)-Y5 

BM(3,15)-Z5 

BM(1,16)-X6 

BM(2,17)-Y6 

BM(3, 18)-Z6 

BM(1,19)-X7 

BM(2,20)-Y7 

BM(3,21)-Z7 

BM(1,22)-X8 

BM(2,23)-Y8 

BM(3,24)-Z8 

BM(4,1)-Y1 

BM(5,2)-Z1 

BH(6,3)-X1 

BM(4,4)-Y2 

BM(5, 5)-Z2 

BM(6, 6)-X2 

BM(4,7)-Y3 
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BM(5,8)-Z3 

BM(6,9)~X3 

BM(4,10)-Y4 

BM(5,11)-Z4 

BM(6,12)-X4 

BM(4,13)-« 

BM(S,14)-Z5 

BM(6,15)-X5 

BM(4,16)-Y6 

BM(5,17)-Z6 

BM(6,18)-X6 

BM(4,19)-Y7 

BM(5,20)-Z7 

BM(6,21)-X7 

BM(4,22)-Y8 

BM(5,23)-Z8 

BM(6,24)-X8 

BH(6,1)-Z1 

BM(4,2)-X1 

BM(5,3)-Y1 

BM(6,4)-Z2 

BM(4,5)-X2 

BM(5,6)-Y2 

BM(6,7)-Z3 

BM(4,8)-X3 

BM(5,9)-Y3 

BM(6,10)-Z4 

BM(4, 11)-X4 

BM(5,12)-Y4 

BM{6, 13)-Z5 

BM(4,14)-XS 

BM(5,15)-Y5 

BM(6,16)-Z6 

BM(4,17)-X6 

BM(5,18)-Y6 

BM(6, 19)-Z7 

BM(4,20)-X7 

BM(5,21)-Y7 

BM(6,22)-Z8 

BM(4,23)-XB 

BM(5,24)-Y8 

RETITRN 

END 


SUBROUTINE GETMS (MM, MN, PROP, S, SAFE) 

DIMENSION S (6), PROP (MM, 10) 

SL-S(l) 

IF(SL.GE.O.O) SLL-PROP(MN,7) 
IF(SL.LT.O.O) SLL-PROP(MN,8) 
SL-ABS(SL) 

SLL-ABS (SLL) 

IF(SLL.EQ.O.O) SLL-1.0 
SAFE- (SLL-SL) /SLL 
D-(S(2)+S(3))/2.0 
B-(S(2)-S(3))/2.0 
B-ABS(B) 

H-S(4)/2.0 

H-ABS(H) 

R-B*B+H*H 
R-SQRT (R) 

MAX-D+R 

MIN*R*”D 

IP(MAX.LE.O.O) GO TO 100 
StL-PROP (MN,9) 

SLL-ABS (SLL) 
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IF(SLL.EQ.O.O) SLL-l.O 
SAF- (SLL-MAX) /SLL 
IF(SAF.LT.SAFE) SAFE-SAF 
100 CONTINUE 

XF(HIN.GE.O.O) GO TO 200 
SLL-PROP (MN, 10) 

SLL-ABS (SLL) 
IF(SLL.EQ.O.O) SLL-l.O 
SAF-(SLL-t-MZN)/SLL 
IF(SAF.LT.SAFE) SAFE-SAF 
200 CONTINUE 
RETURN 
END 


SUBROUTINE GETDD <SSS1 , SSS2, VI, V2, DDl, 002) 

DIMENSION SSS1(6,6) , SSS2 (6, 6) ,001 (6, 6) ,002 (6, 6) 

DIMENSION DD(6,6),CB(6,6) 

DO 7700 11-1,6 
DO 7700 JJ-1,6 

CB (II, JJ) -SSS2 (II, JJ) -SSSl (II, JJ) 

DD(II, JJ)-0.0 
001(11, JJ)-0.0 
7700 002 (II, JJ) -0.0 

DD(1,1)-SSS1(1,1) 

00(1,2)— SSS2(1,1) 

DD(1,3)-SSS1(1,5) 

00(1,4)— SSS2 (1,5) 

DD(1,5)-SSS1(1,6) 

00(1,6)— SSS2(1, 6) 

00(2, D-Vl 
DD(2,2)-V2 
DD(3,1)-SSS1(5,1) 

00(3,2)— SSS2(5,1) 

DD(3,3)-SSS1(5,5) 

00(3,4)— SSS2(5, 5) 

DD(3,5)-SSS1(5,6) 

00(3,6)— SSS2(5, 6) 

DO (4, 3) -VI 
DD(4,4)-V2 
DD(5,1)-SSS1(6,1) 

00(5,2)— SSS2(6,1) 

DD(5,3)-SSS1(6,5) 

00(5,4)— SSS2(6, 5) 

0D(5,5)-SSS1(6,6) 

00(5,6)— SSS2(6, 6) 

DD(6,5)-V1 
00(6, 6) -V2 
CALL INV(DD) 

001 (1,1)-DD(1,2) 

DD1(1,5)-D0(1,4) 

DD1(1,6>-DD(1,6) 

DDl(l,2)-DD(l,l)*CB(l,2)+DO(l,3)*CB(5,2)+DD(l,5)*CB(6,2) 

001 (1,3) -DO (1,1) *CB(l,3)+DO(l,3)*CB(5,3)+DD(l,5)*CB(6,3) 

DD1(1,4)-DD(1,1)*CB(1,4)+DD(1,3)*CB(5,4)+DD(1,5)*CB(6,4) 

001 ( 2 , 2 )- 1.0 

001(3, 3)-1.0 

DD1(4,4)-1.0 

DD1(5,1)-DD(3,2) 

DD1(5,5)-DD(3,4) 

DD1(S,6)-DD(3,6) 

001 (5,2)-DD(3,l)*CB(l,2)+DD(3,3)*CB(5,2)+DD(3,5)*CB(6,2) 
001 (5,3)-DD(3,l)*CB(l,3)+DD(3,3)*CB(5,3)+DD(3,5)*CB(6,3) 
001 (5,4)-DD(3,l)*CB(l,4)+DD(3,3)*CB(5,4)+DO{3,5)*CB(6,4) 
ODl(6,l)-OD(5,2) 

DD1(6,5)-DD(5,4) 
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DD1{6,6)- 

DOl (6,2)< 

DDl (.6,3)' 

DDl (6,4)' 

DD2(1,1)' 

DD2(1,S) 

DD2(1,6) 

DD2(1,2) 

DD2(1,3) 

DD2(1,4) 

DD2(2,2) 

DD2(3,3) 

DD2 (4,4) 

DD2(S,1) 

DD2 (5,5) 

DD2(5,6) 

DD2 (5,2) 

DD2(5,3) 

DD2(5,4) 

DD2 (6,1) 

DD2 (6,5) 

DD2 (6,6) 

DD2 (6,2) 

DD2 (6,3) 

DD2(6,4) 

RETURN 

END 


DD(5'.l *CB(l,2)+DD(5,3)*CB(5,2)*DO(5,5)*CB(6,2) 

SS 5 1 *S(1 3UdD(5,3)*CB(5,3)+DD(5,5)*CB(6,3 
Sd(5;1)*S(1:4UdD(5,3)*<»(5,4)+DD(5,5)*C8(6,4) 

DD(2,2) 

'DD(2,4) 

;m(2:i!*CB(l,2).OD(2,3)*C8(S,2>«TO(2,5)*a<«,J) 

>DD (2# 1) *CB (1# 3) '•'DD (2# 3) *CB (5^3) +DD(2^ 5) ^ 

.DD(2,1)*CB(1,4)+D0(2,3)*CB(5,4)+DD(2,5)*CB(6,4) 

• 1.0 

■ 1.0 

■ 1.0 

-DD(4,2) 

-DD(4,4) 

■•DD (4^1) ’*CB (1,2) +DD (4,3) *CB (5, 2) +DD (4,5) *CB (6,2) 
.SS(1:1)*S(1:3UdD(4,3)*CB(5,3)+DD(4,5)*CT(6,3) 
-DD(4,l)*CB(l,4)+DO(4,3)*CB(5,4)+DD(4,5) CB(6,4) 

-DD(6,2) 

-DD(6,4) 

:?S(^,lNcB(l,2).OD(6,3>-CB(5,2).OD(«,S)*a(6,2l 
-DD(6,1)*CB(1,3)+DD(6,3)*CB(5,3)+DD(6,5)*CT(6,3) 
-DD(6,1)*CB(1,4)+DD(6,3)*CB(5,4)+DD(6,5) CB(6,4) 
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